Quantification of muscle tone

ABSTRACT

A method and device for the quantification of muscle tone, particularly the wrist, wherein non-sinusoidal and non ramp trajectories are used to drive the wrist. Equation 1 is utilized determine the stiffness, viscosity and inertial parameters. 
     
       
         τ s ( t )= K   H θ( t )+ B   H θ( t )+ J   T θ( t )+τ off   [Eq. 1] 
       
     
     wherein where τ s  is the total torque, τ off  is the offset torque, K H  and B H  are the angular stiffness and viscosity of the combined flexor and extensor muscle groups that act on the joint, J T  is the combined inertia of the oscillating appendage, θ is the angular displacement of the system, and          θ   .                   and                   θ   ¨                     
     are the velocity and acceleration. In accordance with the method, the trajectory θ(t) is controlled, the torque response τ(t) is measured and the stiffness, viscosity, and inertial parameters are determined using Equation 1. The method and device are particularly suitable for use on patients with spacisity.

The present application claims the benefit of U.S. provisionalapplication No. 60/230,314, filed on Sep. 6, 2000, incorporated hereinby reference.

FIELD OF THE INVENTION

The present invention relates to a method and apparatus for thequantification of muscle tone. More particularly, the present inventionrelates to a method and device that utilize non-sinusoidal perturbationsto quantify muscle tone. The device is, and method are particularlyuseful in quantifying muscle tone in a spastic patient.

BACKGROUND OF THE INVENTION

Individual skeletal muscle cells are mechanically and anatomicallyarranged in parallel. The total force produced by a muscle is equal tothe sum of the forces generated by its constituent cells. In the normalsubject, muscles that comprise the wrist flexors and extensors arenormally relaxed and are usually recruited to generate force andmovement.

Lower motor neuron paralysis occurs when muscles are deprived of theirimmediate nerve supply from the spinal cord. This occurs when a nervebetween the spinal cord and a muscle is cut or when cell bodies of theventral horn are destroyed as in poliomyelitis. Muscles become soft andatrophic, and reflex response to sensory stimuli is lost. Most nervedisorders that affect limb function are due to upper motor neuronparalysis, wherein damage is present somewhere in the corticospinaltract that originates in the brain and travels through the spinal cord.

Spasticity is defined as abnormal involuntary contraction of a muscle orgroup of muscles due to a rate-dependent reflex mechanism. The spindleelicits the reflex response upon deformation. To a certain extent, thesereflexes are normal and important. In normal operation, these reflexesare suppressed to a certain extent to allow flexibility and motion ofjoints. In spasticity however, there is a disruption in the normalbehavior of the stretch reflex that causes muscles, particularly theflexors, to be extremely resistive to passive stretch (i.e. high intone). As a result, motor control is severely impaired and stiffness ortightness of the muscles may interfere with gait, movement, and speech.Spasticity is usually found in people with some sort of upper motorneuron paralysis, such as those with cerebral palsy, traumatic braininjury, spinal cord injury and stroke patients.

Common symptoms of spasticity may include hypertonicity (increasedmuscle tone), clonus (a series of rapid muscle contractions),exaggerated deep tendon reflexes, muscle spasms, scissoring (involuntarycrossing of the legs) and fixed joints. The degree of spasticity variesfrom mild muscle stiffness to severe, painful, and uncontrollable musclespasms. The condition can interfere with rehabilitation in patients withcertain disorders, and often interferes with daily activities.

Many forms of intervention are available to reduce muscle tone inspasticity. Biochemical pharmaceuticals such as Botox (Botulinum ToxinType A), Intrethecal Baclofen, and Zanaflex (tizanidine) may be used asa biochemical form of intervention. (See R. W. Armstrong, P. Steinbok,D. D. Cochrane, et al. Intrathecally administered baclofen for treatmentof children with spasticity of cerebral origin. J Neurosurg;87(3):409-414, Sep 1997; J. V. Basmajian, K. Shankardass, D. Russel:Ketazolam Once Daily for Spasticity: Double-Blind Cross-Over Study.Arch. Phys. Med. Rehabil. Vol. 67, pp556-557, 1986; P. J. Delwaide:Electrophysiological Analysis of the Mode of Action of Muscle Relaxantsin Spasticity. Annals of Neurology, Vol. 17, No. 1, January 1985, pp90-950). These chemicals are used either to destroy nerve endings at theneuromuscular junction, or they are used as blocking agents whichdepress neuromuscular transmission by competing with acetylcholine forreceptors, thus suppressing nerve conduction. In severe cases,microsurgery is also an option, where incisions are made in thebrainstem or anywhere in the stretch reflex pathway. The safest form ofintervention is physical therapy, such as training, stretching exercisesand casting. (J. C. Otis, L. Root, M. A. Kroll: Measurement of PlantarFlexor Spasticity During Treatment with Tone-Reducing Casts. Journal ofPediatric Orthopedics 5:682-686, 1985).

Tone is defined as the degree of resistance to stretch from an externalsource. An assessment of tone is important in evaluating the degree ofspasticity that a patient has. This assessment is imperative for theclinician to decide what form of intervention to take and to whatdegree. Further, continued assessment throughout intervention isimportant to assess the effectiveness of the intervention. For example,if the Botox dosage administered is too low, it may have little or noeffect in reducing a patient's spasticity. Conversely, if the dosage istoo high, then the patient may lose the ability to control his or herlimb, as blocking too many neuromuscular junctions at the muscle sitemay prevent the central nervous system from having any control over themuscle. An assessment of tone before and after intervention is alsoimportant as it can demonstrate the effectiveness of the treatment.

Probably the most widely accepted clinical test for the evaluation oftone in spasticity is the Ashworth scale shown below in Table 1.

TABLE 1 Grade Description 0 No increase in muscle tone. 1 Slightincrease in muscle tone, manifested by a catch and release or by minimalresistance at the end of the range of motion when the affected part(s)is(are) moved in flexion or extension. 2 Slight increase in muscle tone,manifested by a catch fol- lowed by minimal resistance through theremainder of the range of motion but the affected part(s) is(are) easilymoved. 3 More marked increase in muscle tone through most of the rangeof movement, but affected part(s) easily moved. 4 Considerable increasein muscle tone, passive movement difficult. 5 Affected part(s) is(are)rigid in flexion or extension. TABLE 1: Modified Ashworth scale. TheAshworth scale is a popular, and widely accepted method of evaluatingmuscle tone among clinicians. Though the grades are numeric, they arebased on a qualitative “feel” of resistance to passive stretch movementsperformed by the clinician on the patient. (Taken from Arch. Phys. Med.Rehabili. Vol 80, September 1999).

The clinician moves the subject's limbs about the joints and thenassigns a grade based on a “touchy feely” assessment of how muchresistance the clinician feels. One can easily see the problem here.Since the test is a qualitative one, different clinicians may assigndifferent grades to the same test. Even the same clinician's evaluationmay change due to lack of consistency or depending on whether he or sheis optimistic or pessimistic at the time of the test. The clinician mayalso be biased and may, for example, assign a better grade if he or shehas knowledge of interventions being performed on the patient. Evenputting all these issues aside, it is difficult to get an absolutemeasure of tone using the Ashworth scale. As stated in a review article“The quantification of spasticity has been a difficult and challengingproblem, and has been based primarily on highly observer-dependentmeasurements. The lack of effective measurement techniques has beenrestrictive, since quantification is necessary to evaluate various modesof treatment.” (R. T. Katz, W. Rymer: Spastic hypertonia; mechanisms andmeasurement. Archives of Physical Medicine and Rehabilitation. 1989;70:144-145).

In attempt to quantify muscle tone, some have used Electromyography(EMG) information. (See P. J. Delwaide: Electrophysiological Analysis ofthe Mode of Action of Muscle Relaxants in Spasticity. Annals ofNeurology, Vol. 17, No. 1, January 1985, pp 90-95; A. Eisen:Electromyography in Disorders of Muscle Tone. Le Journal Canadien desSciences Neurologiques, Vol. 14, No. 3. August 1987, pp 501-505; W. G.Tatton, P. Bawa, I. C. Bruce, R. G. Lee: Long Loop Reflexes in Monkeys:An Interpretative Base for Human Reflexes. Cerebral Motor Control inMan.: Long Loop Mechanisms. Prog. clin. Neurophysiol., vol 4, Ed. J. E.Desmedt, pp 229-245, Krager Basel, 1978). EMG electrodes measure andamplify actual action potentials sent to the muscles. Thus, they canmonitor the stretch reflex in action and the active force of muscles canbe monitored.

However, the total force of muscles is a combination of both the activeand passive forces. Active tension is due to muscle stimulation andcontraction due to crossbridge cycling. Independent of musclestimulation and crossbridge cycling, muscles also experience passivetension. Like all materials, muscles experience a passive tension whenthey are stretched beyond their resting lengths. This is due to theinherent mechanical properties of connective tissue in the muscles, suchas elastin. The total force of the muscle is the sum of the active andpassive tensions, as shown in FIG. 1. If muscles were de-innervated,then there would be no active force and the total force would just bethe passive force.

EMG electrodes cannot monitor passive tension. Thus, the EMG signal iscorrelated to the active force exhibited by the muscle but not the totalforce. The most common definition of spasticity is “a motor disordercharacterized by a velocity-dependent increase in muscle tonic stretchreflexes (muscle tone) with exaggerated tendon jerks, resulting fromhyperexcitability of the stretch reflex, as one component of the uppermotor neuron syndrome.” (J. W. Lance: Symposium Synopsis: DisorderedMotor Control, R. G. Feldman, R. R. Young, W. P. Koella, Chicago, YearBook Pulishers, 1980, pp. 485-494). However, it has also been proposedthat an increase in tone is largely caused by changes in the intrinsicmechanical properties of the muscle tissue that causes an increase inthe passive stiffness of the muscle. This contribution of muscle tone isindependent of the stretch reflex and is not encapsulated in the EMGsignal. Still further, EMG signals are extremely noisy, particularly ifsurface electrodes are used. EMG signals measured with surfaceelectrodes can also be influenced by hair, oil, or lotions, and deadskin, thereby i yielding erroneous results. Though the amplitude of theEMG signal is correlated with the active force produced by the muscle,it is extremely difficult to find the proper transformation between thetwo quantities. This relationship differs from person to person due tophysiological differences. It is also sensitive to changes in locationof the electrode within the same subject. Thus, the use of an EMGmachine is not a good option for quantifying muscle tone, but, rather,is more useful for monitoring the timing of reflex responses.

Still other tests, such as the pendulum test, measure the range ofmotion and rate of change of motion of the joint in response to a tendonjerk test. Though such tests may give an indication of tone, thetrajectory of a joint does not give a full picture of muscle tone.

Tone describes the relationship of torque produced by the muscles inresponse to an induced trajectory perturbation, or the resultingtrajectory given a torque perturbation. Either way, the relationshipbetween torque and trajectory can be described using a mathematicalmodel where parameters quantify tone in terms of the viscoelasticproperties of the muscles. However, it is difficult to quantify theseparameters. Unlike heart rate or blood pressure, which are inherentlyphysical quantities that can be measured, tone describes a relationshipbetween torque and displacement, as well as the first derivative ofangular displacement or velocity of the joint based on the viscoelasticproperties of the muscles. Also, inertia of the limb has influence onthe torque, which needs to be properly isolated when determining theviscoelastic properties of the muscles acting on a joint.

In the past, simple DC motors have been used to oscillate an appendageabout its joint while measuring the torque response due to the externalperturbation. The response of appendage movements about the joint ofrotation can be explained in terms of elastic and viscous parameters.These parameters are dependent upon the passive mechanical properties ofmuscles and the active stretch reflex response, which has some inherentdelay. A certain amount of torque is required to move the appendage andto move the parts of the apparatus, as every mass has a rotationalinertia. Thus, the torque response measured can be written as equation 1shown below: $\begin{matrix}{{{\tau_{s}(t)} = {{K_{H}\Delta \quad {\theta (t)}} + {B_{H}{\overset{.}{\theta}(t)}} + {J_{T}{\overset{¨}{\theta}(t)}}}}{{\tau_{s}(t)} = {{K_{H}\left( {{\theta (t)} - \theta_{RP}} \right)} + {B_{H}{\overset{.}{\theta}(t)}} + {J_{T}{\overset{¨}{\theta}(t)}}}}{{\tau_{s}(t)} = {{K_{H}{\theta (t)}} + {B_{H}{\overset{.}{\theta}(t)}} + {J_{T}{\overset{¨}{\theta}(t)}} - {K_{H}\theta_{EQ}}}}{{\tau_{s}(t)} = {{K_{H}{\theta (t)}} + {B_{H}{\overset{.}{\theta}(t)}} + {J_{T}{\overset{¨}{\theta}(t)}} + \tau_{off}}}} & \text{[Eq.~~1]}\end{matrix}$

where τ_(s) is the total torque sensed by the transducer. τ_(off) is theoffset torque and tells the angular position the hand prefers to be inrelative to the origin or bias position. K_(H) and B_(H) are the angularstiffness and viscosity of the combined flexor and extensor musclegroups that act on the joint. J_(T) is the combined inertia ofoscillating appendage, J_(H) and the rotating components of theapparatus, J_(A). J_(T)=J_(H)+J_(A). θ is the angular displacement ofthe system.$\overset{.}{\theta}\quad {and}\quad \overset{¨}{\theta}$

are the velocity and acceleration which are the first and secondderivatives of displacement respectively. θ_(RP) is the angular positionof the resting hand.

This second order differential equation has previously been used as amodel by many including Evans at al (C. M. Evans, S. J. Fellows, P. M.H. Rack, H. F. Ross and D. K. W. Walters: Response of the Normal HumanAnkle Joint to Imposed Sinusoidal Movements. J. Physiol. (1983), 344,pp. 483-502), Rack et al. (P. M. H. Rack, H. F. Ross and T. I. H. Brown:Reflex Response during Sinusoidal Movements of Human Limbs. Prog. Clin.Neurophysiol., vol. 4, Ed. J. E. Desmedt, pp 216-228, 1978), Lehmann etal. (J. F. Lehmann, R. Price, B. J. DeLateur, S. Hinderer, C. Traynor:Spasticity: Quantitative Measurements as a Basis for AssessingEffectiveness of Therapeutic Intervrntion. Arch. Phys. Med. Rehabil. Vol70. pp 6-15, 1989), Prince et al. (R. Price, K. F. Bjornson, J. F.Lehmann, J. F. McLaughlin, R. M. Hays: Quantitative Measurement ofSpasticity in Children with Cerebral Palsy. Developemental Medicine andChild Neurology, 33, pp.585-595, 1991), Minders et al. (M. Meinders, R.Price, J. F. Lehmann, K. A. Questad: The Stretch Reflex Response in theNormal and Spastic Ankle: Effect of Ankle Position. Arch. Phys. Med.Rehabili. Vol 77. pp 487-491, 1996). In each case, simple sinusoidaldisplacements of the ankle or finger joint were used at variousfrequencies. The technique of measuring the frequency response tosinusoidal inputs to investigate properties of second orderelectromechanical systems is common and easy to do. The goal is tocompute K_(H) and B_(H) at various frequencies of oscillation. Since theactivation stretch reflex loop in the central nervous system israte-dependent, one would expect that K_(H) and B_(H) vary at differentperturbation speeds. Sinusouds have been chosen because they are easy togenerate by means of a unidirectional rotating wheel and crank. Nofeedback loop is necessary. It has been argued that since sinusoidalmovements are continuous and smooth, they involve no sudden impulsivemovements that might synchronize a large number of sensory receptors inan artificial or unphysiological way.

However, there is a fundamental mathematical problem in isolating theinertia from the elastic stiffness of the muscles when using simplesinusoidal displacements.

In this case, the displacement, θ(t)=A·sin(ωt), where A is the amplitudeof the sinusoid in radians and ω is the frequency of the oscillation inrads/sec. The velocity,${\overset{.}{\theta}(t)} = {{{\theta}/{t}} = {A\quad {\omega \cdot \cos}\quad \left( {\omega \quad t} \right)}}$

and the acceleration${\overset{¨}{\theta}(t)} = {{{^{2}\theta}/{t^{2}}} = {{- A}\quad {\omega^{2} \cdot \sin}\quad {\left( {\omega \quad t} \right).}}}$

Substituting these state variables into equation 1 and rearranging, weget:

τ_(s)−τ_(off) =K _(H) A sin(ωt)+B _(H) Aω cos(ωt)−J _(T) Aω ² sin(ωt)

τ_(s)−τ_(off)=(K _(H) −J _(T)ω²)(A sin(ωt))+B _(H)ω(A cos(ωt))  [Eq. 2]

The dominant waveform of the torque response is a sinusoid of the samefrequency but φ radians out of phase with the displacement wave as shownin equation 3. After data acquisition, A, M and φ are obtained bylooking at the transfer functions of the displacement and torque signalsin the frequency domain to obtain the magnitude and phases at the forcedfrequency ω

τ_(s)−τ_(off) =M sin(ωt+φ)   [Eq. 3]

From the sum formula:

τ_(s)−τ_(off) =M sin(φ)cos(ωt)+M cos(ω)sin(ωt) If M ₁ =M sin(φ) and M ₂=M cos(φ)

Then:

τ _(s)−τ_(off) =M ₁ cos(ωt)+M ₂ sin(ωt)   [Eq. 4]

Substituting equation 4 into equation 2:

M ₁ cos(ωt)+M ₂ sin(ωt)=(K _(s) −J _(T)ω²)(A sin(ωt))+B _(H)ω(A cos(ωt))

Thus, the total torque measured can be isolated as a component in phasewith the displacement wave and a component quadrature with thedisplacement wave.

M ₁ /A=K _(H) −J _(T)ω²  [Eq. 5]

and

M ₂ /A=B _(H)ω  [Eq. 6]

The amplitude of the in phase component contains both elastic andinertial terms. Researchers have argued that a least squares (error) fitof the in phase component and displacement amplitude with frequencysquared yields a linear relationship of slope intercept form; y=mx+b,where y is M₁/A, x is ω², m is the slope and b is the intercept. Theyhave assumed that the slope is the total inertial and this “totalinertia” value has been used to evaluate the stiffness by manipulatingequation 5 to solve for K_(H) at each frequency.

At first glance of equation 5, this assumption may seem correct. Afterall, the total inertia does remain constant. However, the assumption isfalse because the least squares fit yields an intercept, b, that is asingle approximation of K_(H) across all frequencies. Since the goal isto find a pair of varying stiffness and viscosity values at eachfrequency, this analysis is self-contradictory. This analysis forces theK_(H) values to be trendless (uncorrelated) across frequencies whosevariance depends on how well equation 5 fits to a straight line. To putit simply, there are more unknowns than equations. If N sinusoids areused, each at a different forcing frequency, we are looking for Ndifferent stiffness values and one inertia value. Thus, there are N+1unknowns and N number of equations that resemble equation 5. With moreunknowns than equations, the system is indeterminate and it isimpossible to isolate the stiffness values from the inertia value whenonly sinusoids are used.

Another way to show this problem that prevents isolation of thestiffness and inertial values from the in phase component can be shownusing least squares regression analysis. Usually, we can obtain thelinear parameters of the model in equation 1 by performing a leastsquares regression fit on the data. Given any set of data valuesobtained from a given perturbation trial: $T_{s} = {{\begin{matrix}{\tau_{s}\left( t_{1} \right)} \\{\tau_{s}\left( t_{2} \right)} \\\vdots \\{\tau_{s}\left( t_{n} \right)}\end{matrix}}_{n \times 1}\quad {and}}$ $\Psi = {\begin{matrix}1 & {\theta \left( t_{1} \right)} & {\overset{.}{\theta}\left( t_{1} \right)} & {\overset{¨}{\theta}\left( t_{1} \right)} \\1 & {\theta \left( t_{2} \right)} & {\overset{.}{\theta}\left( t_{2} \right)} & {\overset{¨}{\theta}\left( t_{2} \right)} \\\vdots & \vdots & \vdots & \vdots \\1 & {\theta \left( t_{n} \right)} & {\overset{.}{\theta}\left( t_{n} \right)} & {\overset{¨}{\theta}\left( t_{n} \right)}\end{matrix}}_{n \times 4}$

It is possible to obtain the set of four parameters,$\overset{\sim}{N} = {\begin{matrix}\tau_{off} \\K_{H} \\B_{H} \\J_{T}\end{matrix}}_{4 \times 1}$

that describes the model in equation 1 $\begin{matrix}{{\tau_{s}(t)} = {{K_{H}{\theta (t)}} + {B_{H}{\overset{.}{\theta}(t)}} + {J_{T}{\overset{¨}{\theta}(t)}} + \tau_{off}}} & \text{[eq.1]}\end{matrix}$

using the pseudo-inverse which minimizes the sum squared error

SE=(Ô _(s) −ØÑ)^(T)(Ô _(s) −ØÑ)

V _(emf) =K _(e)

Typically, the pseudo-inverse (denoted by the function pinv in equation7), can be used to obtain the vector of unknowns, P, as follows:

Ñ=pinv(Ø)Ô _(s)=(Ø^(T)Ø)⁻¹Ø^(T) Ô _(s)  [eq. 7]

However, this can usually but not always be done because, in the case ofsinusoidal perturbation at a given frequency, the acceleration waveformis always a constant multiple of the displacement waveform by the scalarquantity −Aω² for all time. Thus the displacement values andacceleration values are linearly dependent. Since the first and thirdrows of Ψ are linearly dependent, the product of the matrices Ψ^(T)Ψ issingular or very close to singular. All singular matrices arenoninvertible, thus P cannot be obtained, or the obtained values areinaccurate.

SUMMARY OF THE INVENTION

The present invention provides a method and apparatus that quantifiesmuscle by using non-sinusoidal perturbations. More particularly, thedevice and method move the upper or lower extremities of a patient in anon-sinusoidal trajectory, θ(t), while measuring the torque responseτ_(s)(t) and utilizing Equation 1 to determine the stiffness K_(H),viscosity B_(H) and inertial J_(T) parameters. $\begin{matrix}{{\tau_{s}(t)} = {{K_{H}{\theta (t)}} + {B_{H}{\overset{.}{\theta}(t)}} + {J_{T}{\overset{¨}{\theta}(t)}} + \tau_{off}}} & \text{[Eq.~~1]}\end{matrix}$

In an exemplary embodiment, as set out herein in more detail, the deviceand method of the present invention can be utilized to quantify theforearm muscle tone, the wrist in particular. The device is particularlyuseful for quantifying muscle tone in a spastic patient. However, it isto be understood that the present method and device are not limited toquantification of forearm muscle tone or to the spastic patient. Rather,the device and method are useful on both the upper and lower extremitiesincluding, for example, the ankle. Further, the method and device can beused to measure muscle tone in the non-spastic patient. For example, themethod and device are useful in evaluating patients with various typesof upper or lower motor neuron paralysis that affects skeletal muscles,such as patients with cerebral palsy, traumatic brain injury, spinalcord injury, stroke, or Parkinson's Disease patients. Further, themethod and device could be used in analyzing a patient's muscle tone inline with, for example, physical therapy that the patient is undergoingto regain control of the muscles in the legs after a spinal accident.

More specifically, in accordance with an exemplary method of the presentinvention, the wrist is driven with an arbitrary trajectory that isneither sinusoidal nor ramp. By controlling the trajectory θ(t) andmeasuring the torque response τ(t), the stiffness, viscosity, andinertial parameters in equation 1 can then be determined.

In general, the device in accordance with the present invention is anautomated tone assessment device that non-invasively and properlyquantifies tone. The device properly determines the muscle tone of aperson's flexors and extensors about a appendage non-invasively. Theappendage is perturbed with arbitrary trajectories by virtue of arobotic design that uses a closed loop automated feedback direct drivedevice that tracks any arbitrary desired trajectory. The device perturbsan appendage with any desired trajectory to properly determine thestiffness of the flexors and extensors of an appendage. Moreparticularly, a desired displacement of the wrist is first determinedusing a set of conditions detailed herein. A discrete set of finitepoints defining this trajectory is then input into a controller, and amotor shaft tracks these points at a certain sampling rate.

Other aspects and embodiments of the invention are discussed infra.

BRIEF DESCRIPTION OF THE DRAWINGS

It should be understood that the drawings are provided for the purposeof illustration only and are not intended to define the limits of theinvention. The foregoing and other objects and advantages of theembodiments described herein will become apparent with reference to thefollowing detailed description when taken in conjunction with theaccompanying drawings in which:

FIG. 1 shows the total muscle tension as the sum of active and passivetensions. (Taken from E. R. Kandel, J. H. Schwartz, and T. M. Jessell.Principles of Neural Science, third edition, Elsevier, New York, 1981,p. 552)

FIG. 2 shows the generation of a smooth filtered random trajectory bylimiting the maximum frequency component of a uniformly distributedrandom number set of 1000 points in the frequency domain. The signal isthen taken back to the time domain via the inverse Fourier transform.The filtered random trajectory is selected to start and stop at theorigin and finally scaled to have an RMS value of 7000 counts (7degrees).

FIG. 3 shows the generation of faster trajectories from the same randomnumber set by increasing the maximum frequency component in thefrequency domain in filtering a random number set.

FIG. 4 shows the generation of two sets of filtered random trajectoriesfrom two sets of random numbers. Trials 1 and 2 have a maximum frequencycomponent of 4 Hz, trials 3 and 4 have a maximum frequency component of5 Hz, trials 5 and 6 have a maximum frequency component of 6 Hz, andtrials 7 and 8 have a maximum frequency component of 8 Hz. Thedisplacements shown are relative to the bias position or starting angleof the wrist. If the left hand is tested instead of the right, thesetrajectories are flipped or negated.

FIG. 5 shows one embodiment of the device in accordance with the presentinvention.

FIG. 6 shows an optical encoder in the form of a code wheel with aseries of slits followed by opaque radial lines. The twoLED/Phototransistor pairs measure the amount of rotation by keepingcount of the number of slits that go by.

FIG. 7 shows the typical circuitry of an LED/phototransistor pair. Whenlight is transmitted from the LED to the phototransistor V_(OUT) isapproximately 5V. When a opaque object blocks the light, V_(OUT) dropsto 0V.

FIG. 8 shows the square wave outputs of two LED/phototransistorchannels. Two channels together are able to detect the direction ofrotation of the code wheel. Two channels also increase the degree ofresolution by 4.

FIG. 9 shows the theory of operation of a simple brushed motor. (Takenfrom M. B Histand, and D. G. Alciatore, Introduction to Mechatronics andMeasurement Systems, WCB/McGraw-Hill, Boston, 1999, p. 315).

FIG. 10 shows the torque-speed curve of the brushless servomotor(B-202-B-31) used in one embodiment of the present invention.

FIG. 11 shows how a patient's arm is properly placed in one embodimentof the device of the present invention.

FIG. 12 shows the proportional relationship between the input voltageand the output current of the amplifier. The constant ofproportionality, or amplifier gain, was experimentally adjusted so thatthe maximum range of voltage (+/−10 V) would correspond to the maximumpossible current output as listed in the manufacturer specificationsheet (+/−6 Amps).

FIG. 13 shows the proportional relationship between the input currentand the torque output of the servomotor. The constant ofproportionality, or torque constant was experimentally determined.

FIG. 14 shows the proportional relationship between the input voltage tothe amplifier and the torque output of the servomotor. The constant ofproportionality between the two is the product of the torque constantand the amplifier gain.

FIG. 15 is a block diagram showing components of the PD closed loopcontrol system used to follow the desired trajectories in accordancewith one embodiment of the present invention.

FIG. 16 shows a typical analytical frequency response and step responseof the entire closed loop system shown in FIG. 15.

FIG. 17 shows four possibilities of bias positions that subjects weretested at in accordance with the present invention.

FIG. 18 shows raw data collected from optical encoder, torque transducerand EMG electrodes in trial 1 of experiment DMJF.

FIG. 19 shows raw data collected from optical encoder, torque transducerand EMG electrodes in trial 2 of experiment DMJF.

FIG. 20 shows raw data collected from optical encoder, torque transducerand EMG electrodes in trial 3 of experiment DMJF.

FIG. 21 shows raw data collected from optical encoder, torque transducerand EMG electrodes in trial 4 of experiment DMJF.

FIG. 22 shows raw data collected from optical encoder, torque transducerand EMG electrodes in trial 5 of experiment DMJF.

FIG. 23 shows raw data collected from optical encoder, torque transducerand EMG electrodes in trial 6 of experiment DMJF.

FIG. 24 shows raw data collected from optical encoder, torque transducerand EMG electrodes in trial 7 of experiment DMJF.

FIG. 25 shows raw data collected from optical encoder, torque transducerand EMG electrodes in trial 8 of experiment DMJF.

FIG. 26 shows raw data collected from optical encoder, torque transducerand EMG electrodes in trial 9 of experiment DMJF.

FIG. 27 shows raw data collected from optical encoder, torque transducerand EMG electrodes in trial 10 of experiment DMJF.

FIG. 28 shows the frequency response of the Savitsky-Golay filter usedon the encoder counts.

FIG. 29 shows the frequency response of the butterworth filter used onthe torque signal.

FIG. 30 shows rectified data filtered from the optical encoder, torquetransducer and EMG electrodes in trial 1 of experiment DMJF.

FIG. 31 shows rectified data filtered from the optical encoder, torquetransducer and EMG electrodes in trial 2 of experiment DMJF.

FIG. 32 shows rectified data filtered from the optical encoder, torquetransducer and EMG electrodes in trial 3 of experiment DMJF.

FIG. 33 shows rectified data filtered from the optical encoder, torquetransducer and EMG electrodes in trial 4 of experiment DMJF.

FIG. 34 shows rectified data filtered from the optical encoder, torquetransducer and EMG electrodes in trial 5 of experiment DMJF.

FIG. 35 shows rectified data filtered from the optical encoder, torquetransducer and EMG electrodes in trial 6 of experiment DMJF.

FIG. 36 shows rectified data filtered from the optical encoder, torquetransducer and EMG electrodes in trial 7 of experiment DMJF.

FIG. 37 shows rectified data filtered from the optical encoder, torquetransducer and EMG electrodes in trial 8 of experiment DMJF.

FIG. 38 shows rectified data filtered from the optical encoder, torquetransducer and EMG electrodes in trial 9 of experiment DMJF.

FIG. 39 shows rectified data filtered from the optical encoder, torquetransducer and EMG electrodes in trial 10 of experiment DMJF.

FIG. 40 shows the amount of residual torque when using one set of timeinvariant parameters, {tilde over (τ)}_(off), {tilde over (K)}_(H),{tilde over (B)}_(H), and {tilde over (J)}_(T) over the entire data setin each trial in experiment DMJF.

FIG. 41 shows the results for control adult subject JF. JF's right armwas tested at a bias of both 0 degrees and 30 degrees flexion at each ofthe ten filtered, random trajectories.

FIG. 42 shows results for control adult subjects PW, JM, and DD. One armwas tested at a bias of both 0 degrees and 30 degrees flexion at each ofthe ten filtered, random trajectories.

FIG. 43 shows results for control adult subjects JF, MT, and RP. One armwas tested at a bias of both 0 degrees and 30 degrees flexion at each ofthe ten filtered, random trajectories.

FIG. 44 shows results for control adult subjects SN, JC, and SB. One armwas tested at a bias of both 0 degrees and 30 degrees flexion at each ofthe ten filtered, random trajectories.

FIG. 45 shows results for control adult subjects EK, SW, and JH. One armwas tested at a bias of both 0 degrees and 30 degrees flexion at each ofthe ten filtered, random trajectories.

FIG. 46 shows results for control adult subjects SP, VG, and NM. One armwas tested at a bias of both 0 degrees and 30 degrees flexion at each ofthe ten filtered, random trajectories

FIG. 47 shows results for control adult subjects SL, SR, and HH. One armwas tested at a bias of both 0 degrees and 30 degrees flexion at each ofthe ten filtered, random trajectories.

FIG. 48 shows results for control adult subjects LA, CS, and SD. One armwas .tested at a bias of both 0 degrees and 30 degrees flexion at eachof the ten filtered, random trajectories.

FIG. 49 shows results for control adult subjects RB, LA, and HC. One armwas tested at a bias of both 0 degrees and 30 degrees flexion at each ofthe ten filtered, random trajectories.

FIG. 50 shows results for control adult subjects KT, WZ, and MH. One armwas tested at a bias of both 0 degrees and 30 degrees flexion at each ofthe ten filtered, random trajectories.

FIG. 51 shows results for control adult subjects LH, DR, and RW. One armwas tested at a bias of both 0 degrees and 30 degrees flexion at each ofthe ten filtered, random trajectories.

FIG. 52 shows results for control adult subject GD. One arm was testedat a bias of both 0 degrees and 30 degrees flexion at each of the tenfiltered, random trajectories.

FIG. 53 shows the average {tilde over (τ)}_(off), {tilde over (K)}_(H)and {tilde over (B)}_(H) values across all control adults for eachtrail. The eleventh bar shows the overall average and standard deviationacross all control adults, across all ten trials.

FIG. 54 shows the results for adult subjects with hemiplegic spasticity:AM, JS, and SH. The affected and unaffected arms were tested at a biasof 30 degrees flexion at each of the ten filtered, random trajectories.

FIG. 55 shows the results for adult subject BY with hemiplegicspasticity. The affected and unaffected arms were tested at a bias of 30degrees flexion at each of the ten filtered, random trajectories.

FIG. 56 shows the average {tilde over (τ)}_(off), {tilde over (K)}_(H)and {tilde over (B)}_(H) values across all adults with hemiplegicspasticity for each trail. The eleventh bar shows the overall averageand standard deviation across all adults with spasticity, across all tentrials.

FIGS. 57-124 show block diagrams of code used in accordance with oneembodiment of the present invention.

FIG. 125 shows one embodiment of the wiring schematic of the motorsystem in accordance with the present invention

FIG. 126 shows one embodiment of the wiring schematic from the sensorsto the data acquisition board in accordance with the present invention.

FIG. 127 shows the pulse width of the embedded chip in the US Digitaldecoder was decreased via the addition of a resistor as shown. This wasdone because the original pulse width was wide and the optical encoderwas skipping counts at high rotational speeds.

FIG. 128 shows an AutoCAD mechanical design drawing of one embodiment ofthe motor coupling in accordance with the present invention.

FIG. 129 shows an AutoCAD mechanical design drawing of one embodiment ofthe shaft coupling in accordance with the present invention.

FIG. 130 shows an AutoCAD mechanical design drawing of one embodiment ofthe output shaft in accordance with the present invention.

FIG. 131 shows an AutoCAD mechanical design drawing of one embodiment ofthe forked wrench in accordance with the present invention.

FIG. 132 shows an AutoCAD mechanical design drawing of one embodiment ofthe hand support in accordance with the present invention.

FIG. 133 shows an AutoCAD mechanical design drawing of one embodiment ofthe hand plate in accordance with the present invention.

FIG. 134 shows an AutoCAD mechanical design drawing of one embodiment ofthe poles that secure the hand on the hand plate in accordance with thepresent invention.

FIG. 135 shows an AutoCAD mechanical design drawing of one embodiment ofthe arm rest attachment in accordance with the present invention.

FIG. 136 shows an AutoCAD mechanical design drawing of one embodiment ofthe Mechanical stop pin in accordance with the present invention.

FIG. 137 shows an AutoCAD mechanical design drawing of one embodiment ofthe top surface of the housing in accordance with the present invention.

FIG. 138 shows an AutoCAD mechanical design drawing of one embodiment ofthe bottom surface of the housing in accordance with the presentinvention.

FIG. 139 shows an AutoCAD mechanical design drawing of one embodiment ofthe back surface of the housing in accordance with the presentinvention.

FIG. 140 shows an AutoCAD mechanical design drawing of one embodiment ofthe left side surface of the housing in accordance with the presentinvention.

FIG. 141 shows an AutoCAD mechanical design drawing of one embodiment ofthe right side surface of the housing in accordance with the presentinvention.

FIG. 142 shows an AutoCAD mechanical design drawing of one embodiment ofthe encoder mount in accordance with the present invention.

FIG. 143 shows an AutoCAD mechanical design drawing of one embodiment ofthe motor mounting plate in accordance with the present invention.

FIG. 144 shows an AutoCAD mechanical design drawing of one embodiment ofthe front surface of the housing in accordance with the presentinvention.

FIG. 145 shows an AutoCAD mechanical design drawing of one embodiment ofthe bearing on the top surface of the housing in accordance with thepresent invention.

FIG. 146 shows an AutoCAD mechanical design drawing of one embodiment ofthe motor support column in accordance with the present invention.

FIG. 147 shows an AutoCAD mechanical design drawing of one embodiment ofthe armrest slide in accordance with the present invention.

FIG. 148 shows an AutoCAD mechanical design drawing of one embodiment ofthe bearing clamp ring in accordance with the present invention.

FIG. 149 shows an AutoCAD mechanical design drawing of one embodiment ofthe hand sub assembly in accordance with the present invention.

FIG. 150 shows an AutoCAD mechanical design drawing of one embodiment ofthe complete assembly drawing in accordance with the present invention.

FIG. 151 shows an AutoCAD mechanical design drawing of one embodiment ofthe labeled parts in accordance with the present invention.

DETAILED DESCRIPTION OF THE INVENTION

The present invention provides a method for quantifying muscle tone andan apparatus that carries out the method.

In general, the device in accordance with the present invention is anautomated tone assessment device that non-invasively and properlyquantifies tone. Quantification of tone allows a clinician to decide onthe nature and degree of intervention to administer. A pre andpost-evaluation can also determine the effectiveness of the interventionas well as track the progress of the patient. The system software caneasily be modified to output the muscle stiffness of the patient to theclinician in real time. The implications of this are in the operatingroom where a neuro-surgeuon makes temporary lesions in the brain whilemonitoring the patient's muscle tone in real time. Once the surgeonfinds the proper location in the brain that alleviates the degree ofmuscle tone with out too much loss in function, he/she can make thelesion permanent.

In an exemplary embodiment the device and method of the presentinvention are utilized to quantify the forearm muscle tone, inparticular, the wrist in a spastic patient. However, it is to beunderstood that the present method and device are not limited toquantification of forearm muscle tone or to the spastic patient.

When used to quantify muscle tone in the spastic patient, upondetermination of the muscle tone, the degree of spasticity that thepatient has can be determined, thereby allowing a clinician to decidewhat form of intervention to take and to what degree. The muscle tonecan further be monitored over time to allow the clinician to determinewhether the intervention is providing benefits to the patient.

More particularly, the present invention utilizes equation 1 todetermine the stiffness, viscosity and inertial parameters.$\begin{matrix}{{\tau_{s}(t)} = {{K_{H}{\theta (t)}} + {B_{H}{\overset{.}{\theta}(t)}} + {J_{T}{\overset{¨}{\theta}(t)}} + \tau_{off}}} & \text{[Eq.~~1]}\end{matrix}$

While equation 1 has been used in the past, the methods and devices thatwere used to apply equation 1 utilized sinusoidal displacements of theankle or finger joint at various frequencies. Using sinusoidaldisplacements, as set out above, leads to an indeterminate system, andit is impossible to isolate the stiffness values from the inertiavalues.

In accordance with the method of the present invention, a desiredtrajectory that is non-sinusoidal and not a ramp is first determined.The desired trajectory is then input into the device of the presentinvention, which is an electromechanical system that utilizes a feedbackcontroller to track the desired displacement trajectory.

In accordance with the present method, a desired displacement of thewrist is first determined using the following set of conditions.

(1) The trajectory is ‘smooth’. In other words, when perturbing anappendage about its joint of rotation, there should be “no suddenimpulsive movements that might synchronize a large number of sensoryreceptors in an artificial or unphysiological way.” (See P. M. H. Rack,H. F. Ross and T. I. H. Brown: Reflex Response during SinusoidalMovements of Human Limbs. Prog. Clin. Neurophysiol., vol. 4, Ed. J. E.Desmedt, pp 216-228, 1978). As used herein, a “smooth discretetrajectory” is mathematically defined as follows:

θ_(d)(t_(k)) is a discrete signal consisting of set of N (1000) pointsthat could be sampled from a continuous, differentiable function,θ_(d)(t) for all time t such that:

a. The sampling frequency (fhz_(mc)) of the discrete signal must be nearto or greater than or equal to 100 Hz, particularly if the maximumfrequency required to represent the discrete signal approaches orexceeds 10 Hz.

b. In the frequency domain, the maximum frequency needed to representthe desired trajectory must be less than the smaller of one-sixthfrequency of the sampling rate of the motor controller, fhz_(mc), or 15Hz.

(2) At time, t=0, θ_(d)(0)=0, the motor does not jump rapidly to theinitial point.

(3) The trajectory θ_(d)(t) is periodic (repeatable) and smooth for alltime t>0. In other words θ_(d)(t_(N+1)) equals θ_(d)(t₁).

(4) The trajectory has a finite range of motion with a controllable,defined distribution.

(5) The matrix Ψ has linearly independent columns, which ensures thatthe matrix Ψ^(T)Ψ is nonsingular.

Filtered random trajectories are then generated using, for example,MATLAB in the following way. There exist three discrete signals,r_(a)(t_(k)), r_(b)(t_(k)) and r_(c)(t_(k)) in the time domain whichcorrespond to R_(a)(ω), R_(b)(ω), and R_(c)(ω) respectively in thefrequency domain. The time between two consecutive points given to thecontroller is dt_(mc) seconds apart.

(1) A set of N(1000) uniformly distributed numbers is generated from−0.5 to 0.5. The frequency of the motor controller, fhz_(mc), is 128 Hz,so consecutive points are dt_(mc)=0.0078 seconds apart in time. Thisdiscrete signal is called r_(a)(t_(k)).

(2) The discrete Fourier Transform of r_(a)(t_(k)) is taken, R_(a)(ω) toobtain the magnitude, |R_(a)(ω)| and phase, φ(R_(a)(ω)) of the signal inthe frequency domain.

(3) Now, R_(b)(ω) is generated. For all ω, φ(R_(b)(ω))=φ(R_(a)(ω)). Anω_(max) is selected such that |R_(b)(0)|=0, for 0<ω≦ω_(max),|R_(b)(ω)|=|R_(a)(ω)|, and for all ω>ω_(max), |R_(b)(ω)|=0. In FIG. 2,ω_(max)=2π*4 Hz.

(4) The Inverse Fourier Transform of R_(b)(ω) is taken to obtainr_(b)(t). Note that the mean value of r_(b)(t_(k))=0.

(5) There is no guarantee that r_(b)(0₁)=0, however, the trajectory doescross θ=0 at several instances in time. The index, min_i is found wherer_(b)(tmin_i) has the lowest absolute value which is approximately 0.Note that r_(b)(t_(k)) is a continuous periodic function for all timewith a period of T_(θ)=dt_(mc)*(N)=7.8125 s.

(6) r_(c)(t_(k))=α(r_(b)(t_(k)+t_(min) _(—) _(i))−r_(b)(t_(min) _(—)_(i))). By shifting the signal r_(b)(t_(k)) back t_(min) _(—) _(i)seconds in time, r_(c)(0)=0. α is a scalar number which scales the rangeof motion of the trajectory. α was selected so that the root meansquared value of r_(c)(t_(k)) is 7000/{square root over (2)} counts.Thus, α=7000*RMS(r_(b)(t_(k)))/{square root over (2)}

(7) Finally, each data point in r_(c)(t_(k)) is rounded to the nearestinteger which gives the desired trajectory, θ_(d)(t_(k)). This is donebecause the encoder (position sensor) and the controller deal with aninteger number of digital counts.

The discrete set of finite data points defining the trajectory is theninputted into the controller, and a motor shaft tracks these points at acertain sampling rate. The desired trajectories generated in thisfashion satisfy all conditions listed above. Faster desired trajectoriescan be further generated that cover the same distribution of angularpositions from the same random number set. For example, as shown in FIG.4, two random number sets were used with maximum frequency components of4, 5, 6, 7, and 8 Hz to give 10 trials.

Generally, the trajectory is determined for testing of the right hand.However, if the left hand is to be tested, these trajectories for theright hand can simply be flipped or negated about the bias position, orrelative origin, 0, to ensure symmetry across contralateral limbs.

The device in accordance with the present invention is anelectromechanical system that is designed to drive the wrist with thefiltered random trajectories described above. Because these trajectorieshave complex variable position and velocity profiles, an automatedfeedback controller, or robot, is used which constantly varies thetorque, speed and direction of the actuator (motor). The desiredtrajectory is given as input, and the actuator moves back and fourthtracking the desired trajectory. Thus, design of the system inaccordance with the present invention includes a mechanical design, theuse of sensors and actuators, feedback control, data acquisition andprogramming.

One embodiment of the electromechanical system is shown in FIG. 5. Asshown, the device includes a housing 1. The materials used infabricating the housing are not particularly limited provided theysupply the required structural support for holding a patient's arm onthe top surface 2 of the housing 1 during testing. One preferredmaterial is a metal such as aluminum. The housing 1 is shown as box-likein shape, however, the shape of the housing 1 is not limited and cancome in a variety of shapes. Preferably, the top surface 2 of thehousing 1, where the patient's arm rests during testing, is flat.

On the top surface 2 of the housing 1 is a hand securing mechanism 3where the patent places his hand during the test. The hand securingmechanism 3 is designed such that the patient's hand can be placedsideways within the hand securing mechanism 3 with the palm facing theleft (when the right hand is tested) and the thumb facing upwards. Toaccommodate various sized hands, the hand securing mechanism 3 ispreferably adjustable to varying sizes.

In one preferred embodiment, as shown in FIG. 5, the hand securingmechanism 3 is formed of two vertical bars 6 that extend upwards fromthe top surface 2 of the housing 1 so that the patient can place his orher hand between the two bars as shown in FIG. 7. The two vertical bars6 preferably have padding or cushions 7 on their outer surfaces toprovide comfort, as shown in FIG. 11. The two vertical bars 6 mayaccommodate various sized hands by, for example, slidably attaching thetwo vertical bars 6 to the housing 1 so that they can be moved inwardsand outwards and locked into a desired position. The padding or cushions7 can also be removable and can come in a variety of thicknesses toprovide smaller and larger openings between the two vertical bars 6,thereby accommodating various sized hands. In yet another embodiment,the padding or cushions 7 are compressible and expandable such that thepadding or cushions 7 compress or expand as required to hold hands ofvarious sizes in place.

A hand rest 4 and an arm rest 5 may further be mounted on top surface 2of the housing 1 to assist in providing proper positioning of the armand hand during testing. Preferably, the arm rest 5 includes an armsecuring mechanism 10 such as, for example, one or more Velcro straps,as shown in FIG. 11, to secure the arm in place during testing.

In a preferred embodiment, a motor shaft 11 drives the hand about thewrist joint back and forth to follow the desired trajectories. In thisembodiment, a direct drive system is preferably used wherein the handrest 4 is directly connected to an extension of the motor shaft 11,called the output shaft 12. Preferably, there are no gears in betweenthe motor shaft 11 and the output shaft 12 so that there is nojitteriness or backlash when the motor shaft 11 moves or changesdirection. During use, the hand is positioned on the hand rest 4, whichmay be in the form of a small metal platform as shown in FIG. 6. Thehand is secured as it rotates about the wrist joint.

The hand securing mechanism 3, e.g. two vertical bars 6, secure the areaof the hand directly proximal to the first joint between the metacarpalsand phalanges so that the knuckles lie past the hand securing mechanism3, as shown in FIG. 11. The hand is secured so that the patient beingtested does not grasp the hand securing mechanism 3. In fact, thepatient should be completely relaxed so that all muscle contractions arecompletely involuntary due to the stretch reflex arc. Thus, the patientshould not grasp the hand securing mechanism 3 during the test, as thereare muscle groups that control flexion/extension of the wrist as well asthe fingers.

In a preferred embodiment, the hand rest 4 lies on a sliding forkedwrench 22 that is screwed into the output shaft 12. The hand rest 4 isthen secured to a position on a sliding track, such that the axis ofrotation of the wrist joint lies directly above the axis of rotation ofthe output shaft 12. During the perturbation of the hand, the remainderof the arm is secured to the arm rest 5 with the arm securing mechanism10 as set out above. The arm rest 5 is preferably made out of lowtemperature thermoplastic that is covered with a thin layer of foam forcomfort. In a preferred embodiment, the arm rest 5 is on an adjustabletrack, e.g. a friction screw track, that allows for minor adjustments inthe positioning of the arm and hand if necessary.

The arm rest 5 and associated arm securing mechanism 10 hold the armsecurely in place with quick setup time and is sized to accommodatevarious arm sizes. Removable cushions or pads of varying sizes can, forexample, be placed in-between the arm rest 4 and the arm to accommodatevarious sized arms. In one embodiment, a compressible and expandablecushion (not shown) is mounted on the arm rest 5 such that the cushionfits about varying sized arms by compressing or expanding as necessary.Further, the arm rest 5 and arm securing mechanism 10 is preferablyadjustable to varying sized arms. For example, in one embodiment, Velcrostraps, which can be quickly and easily secured tighter or looser, areused as the arm rest arm securing mechanism 10. In another embodiment,one or more straps having multiple snaps along the length of the strapscan be used for varying sized fits. Still further, one or more strapssimilar to a belt with multiple holes can be used to provide varyingsized fits.

A torque transducer 20 (sensor) is placed between the motor shaft 11 andthe output shaft 12. The torque transducer 20 may, for example, beplaced between the motor shaft 11 and the output shaft via two machinedcouplings secured with screws. Preferably, a top coupling has a flatsurface for a flat end screw to securely hold the output shaft 12 ontothe coupling. During fabrication of the device, before the top surface 2of the housing 1 is put in place, an incremental optical encoder 24(position sensor) is placed around the output shaft 12 as shown in FIG.5. After the top surface 2 is secured on top of the housing 1, the innerrotating ring 26 of the encoder 24 is secured to the output shaft 12with, for example, a clamp adapter ring positioned between the opticalencoder 24 and output shaft 12. In a preferred embodiment, the rotatingring 26 is secured to the output shaft 12 via a plurality of hex-screwswith the clamp adapter ring in between the inner rotating ring 26 of theoptical encoder 24 and the output shaft 12. A rotation of the outputshaft 12 corresponds to the same amount of rotation of the code wheeldisc of the encoder 24.

The incremental optical encoder 24 is a device that converts motion intoa series of digital pulses that can be measured. A preferred incrementaloptical encoder 24 in accordance with the present invention is shown inFIG. 6 and consists of a disc with a series of equally spaced slits ortransparent segments 32 on it. Separating these slits or segments 32 areopaque radial lines 38 of equal arc thickness. On either side of thedisk are two LED emitters 34 paired with two phototransistors 36. Thetypical circuitry of an LED/phototransistor pair is shown in FIG. 7.

Typically the forward voltage drop across the LED, V_(AK) is rated ataround 1V at 20 mA. Since V_(S) is 5V, the resistor, R₁ is needed todissipate 4V. Thus by Ohm's law,

R ₁=4V/20 mA=200Ω

When there is no opaque object between the emitter and detector, thelight from the emitter supplies a base current to the transistor. ByKirchoff's Voltage Law,

V_(C)=V_(S)=5V

The saturated collector emitter voltage, V_(CE) is rated at 0.2V, thus

V_(out)=5V−0.2V=4.8V

If R₂=1MΩ,

I _(E)=V_(out) /R ₂=4.8V/1MΩ=4.8 μA

By Kirchoff's Current Law, inside the transistor the emitter current isequal to the sum of the base and collector currents,

I _(E) =I _(C) +I _(B)

Also in the transistor,

I _(C) =h*I _(B)

where h is some positive number. In summary, when there is no opaqueobject between the LED emitter 34 and the phototransistor 36, there issome finite base current, I_(B) from the LED emitter 34, and V_(out) hasa reading of almost 5V (TTL high). When the opaque radial line blocksthe light between the LED emitter 34 and the phototransistor 36,I_(B)≈0, and I_(C)≈0, thus V_(out)≈0V. As the disc 30 rotates betweenthe transparent slits 32 and the radial lines 38, the output voltageshifts between TTL high and TTL low respectively, resulting in a squarewave output. Typically, the number of square wave cycles (ncypr)generated by an LED/phototransistor pair 34, 36 is equal to the numberof radial lines 38 or transparent slits 32 on the disc 30. However ifmore sophisticated circuitry and interpolation techniques is used, thenumber of square wave cycles may be greater than the number of radiallines 38 on the disc. The preferred incremental optical encoder used inthe present device is the Gurley hollow shaft encoder, which has 11,250radial lines on the disc and generates 90,000 square wave cycles perrevolution.

The more radial lines 38 (square wave cycles) that there are on the disc30, the higher the count resolution. The simplest way to obtain angularposition from a pair of LED emitters 34 and phototransistors 36 is tocount the number of cycles that go by, either by counting the rising orfalling edges of one square wave output (CH A) as depicted in FIG. 8,(c) and (d). In this case, the angular resolution is equal to the numberof cycles per revolution [ncypr] counts/rev or [360/ncypr] degrees. Asshown in FIG. 8, (e) and (f), both the rising and falling edges of CHAare counted, the resolution can be doubled to [2*ncypr] counts/rev or[720/ncypr] degrees.

If only one channel of square wave output is used, it is impossible totell the direction of rotation of the disc 30 as the counting pattern ofthe clockwise and counterclockwise direction are the same (see FIGS.9(c) and 9(d)). To overcome this problem, an additionalLED/phototransistor pair 34, 36 is physically placed around the disk 30such that it's square wave output (CH B) is a quarter cycle out of phasewith that of CH A. Thus, the square wave cycles from both channels arecalled quadrature signals. When the disk 30 rotates in one direction(i.e. clockwise), CH A leads CH B by a quarter of a cycle, but when thedisk 30 rotates in the opposite direction, CH A lags CH B by a quarterof a cycle. Another advantage of using quadrature signals is that thecount resolution of the optical encoder 24 can be increased to [4*ncypr]counts/rev by counting the rising and falling edges of both CH A and CHB. This is done using a quadrature decoder chip (for example, theUS-Digital PC-6-84-4). As shown in FIG. 8, (g) and 5(h), the countpatterns in the clockwise and counterclockwise directions aredifferent—they are reversed. The quadrature decoder converts the signalsfrom CH A and CH B into a direction (up/down pulse) and count pulses,where a count pulse is generated every time the TTL state of eitherchannel changes. Since the Gurley hollow shaft encoder, preferably usedin the present device, generates 90,000 cycles per revolution in eachchannel, after quadrature edge detection, the encoder gives 360,000counts per revolution, which is an amazing resolution of one thousandth(10⁻³) of a degree.

In the setup of the present invention, the outputs of the opticalencoder 24 split and go to both the motion controller for feedbackcontrol and to a National Instruments data acquisition card foracquiring data. Theoretically, position information can be obtained byinterrogating the motion controller via software, in which case, dataacquisition of the position signal is not necessary. However, when themotion controller is interrogated at a high sampling rate, its operationand control of the motor shaft 11 is impeded. The controller softwarecannot both operate the feedback loop to follow a desired trajectory andoutput data to the user at 500 Hz.

The optical encoder 24 outputs three additional signals relative toground, not just A and B. These outputs are the index signals,{overscore (A)}, and {overscore (B)}. Once activated, the opticalencoder 24 keeps track of angular position, in counts, relative to theinitial angular position where it started counting. The index signalprovides an absolute reference so that the absolute position of the handrest 4 is known in space. In addition to the many equally spaced slitson the disk 30, there is one additional slit that is at a differentradial distance away from its center, where the index signal is at TTL 1only in one absolute position using an additional LED/phototransistorpair 34, 36 as described above. The counter value determined by channelsA and B is reset to zero where the TTL value of the index is 1 and,thus, the index serves as an origin in space. Once this position isfound, the system knows where the absolute position of the hand rest 4is in real space.

As used herein, the axis of rotation of the output shaft 12 is definedas pointing upward, so that a counter clockwise rotation is positivechange in direction, and a clockwise rotation is negative change indirection. The motion controller however, defines a clockwise rotationas positive, and a counterclockwise rotation as negative.

Preferably, in the device of the present invention, the index signal isnot used. Rather, in one preferred embodiment, the reference point usedto determine absolute position is the left mechanical stop 26. A rightmechanical stop 27 may further be located opposite the left mechanicalstop 26. The motor shaft 11 moves counterclockwise at a slow speed untilthe forked wrench of the hand rest 4 hits the left mechanical stop 26 atwhich point the motor stops. In this configuration, the mid-axis of theforked wrench is 90 degrees (90940 counts) away from the midline, thecounter value with the motion controller is set to −90940 which is−90.940 degrees. Now the counter value reads, and will continue to read,the correct absolute position until the computer that powers the opticalencoder 24 is turned off.

The output pulse width determines the minimum separation betweenconsecutive A and B pulses that can be detected by the decoder chip. Ifthe separation between consecutive A and B pulses is less than theoutput pulse width of the decoder, then two output pulses overlap andonly one count is registered instead of two, resulting in an erroneousposition reading. The width of the clock pulse is influenced by thevalue of the resistor between pins 1 and 3 of the LS7084 quadraturedecoder chip embedded in the US-Digital PC-6-84-4 as shown in FIG. 127.In an exemplary embodiment, an added resister across pins 1 and 3 isused to decrease the decoder pulse width to less than 750 ns.

An actuator is a device that converts electromechanical force intomotion. This electromechanical force can be explained by Lorentz's Law:

{right arrow over (F)}={right arrow over (I)}×{right arrow over (B)}

When a current in a conductor, {right arrow over (I)}, is exposed to amagnetic field, {right arrow over (B)}, a force {right arrow over (F)}is produced in the direction perpendicular to both {right arrow over(I)} and {right arrow over (B)} as described by the right hand rule.

Since the hand is perturbed about the wrist joint during testing, amotor 40 was used as the actuator for this task. Generally, thestationary outer housing of a motor, the stator, is magnetized either bypermanent magnets or by currents flowing through coiled wires wrappedaround iron cores. The rotor or rotating shaft is also magnetized eitherby permanent magnets or winded field coils. The rotor and its windingsare collectively known as armature. If the rotor is magnetized viacoils, then it has commutator segments where the current is supplied.Brushes provide stationary electrical contact between the power sourceand the moving commutator and rotor. FIG. 9 shows the theory ofoperation of a simple two-pole DC motor in view of the long axis of therotor, wherein “DC” denotes Direct Current. The setup shown consists ofa stator that is a permanent magnet and a rotor that is magnetized viacoils. In the starting position, (i), the right brush touches commutatorsegment A and the left brush touches segment B creating a current thatflows through the rotor winding. This current magnetizes the rotor intoN and S poles as shown. Since like poles repel and opposite polesattract, the rotor rotates in the clockwise direction as shown in (ii)and (iii). In (iv), for a brief instant, the commutators lose contactwith the brushes, no current flows through the armature and the rotor isdemagnetized. However, the rotor continues to turn clockwise due tomomentum. In (v), the commutator contacts switch brushes and thearmature current now flows in the opposite direction, so the poles ofthe rotor are switched. Thus, the motor continues rotating as thearmature current constantly switches direction, maintaining torque inthe clockwise direction. Like most motors, the motor shown in FIG. 9 isa brush motor. Its stator is a permanent magnet and its rotor armaturehas a constantly switching current provided by an external currentsource through the brushes. The brushes of a DC motor have severallimitations: brush life, brush residue, maximum speed, and electricalnoise. The currents produced in the armature produce a large amount ofheat in the rotor that is difficult to dissipate.

The brushless DC motor operates by the same principle as the brush motorexcept that the rotor consists of a permanent magnet and the stator ismagnetized via rotating fields in the stator. It is like a brush motorturned inside out. Indicative of its name, the brushless motoreliminates the need for brushes and commutators. Brushless DC motors arepotentially cleaner, faster, more efficient, less noisy, more reliable,produce no sparks and the rotor does not heat up. These motors are alsoknown as permanent magnet motors. Motors that use feedback for aposition or trajectory control applications are called servomotors.

Preferably, the motor according to the present invention is a brushlessDC servomotor, although other types of motors are compatible for use inthe present invention. A particularly preferred motor for use in thepresent invention is the B-202-B-3 1 brushless servomotor fromKollmorgen as its speed range, torque-speed curves and maximum stalltorque are acceptable for use in the present invention. The speed torquecurve shows the maximum torque the motor can produce at a certain speedusing the amplifier of the present invention (see FIG. 10). For thepresent invention, this is the torque applied by the hand in response tothe perturbations given, and, more importantly, the inertias of the handand the hand rest. As shown in FIG. 10, the maximum possible operatingspeed of the motor is 3800 RPM. The motor can generate torques up to 3.5Nm while maintaining a speed fairly close to the maximum value. If themotor is to generate torques above 3.5 N, the maximum speed it canmaintain drops very sharply until the motor finally stops moving at themaximum torque of 4.7 Nm. This maximum torque value where the motorshaft is stationary is called the stall torque.

When the load torque on the motor is equal and opposite to the toqueproduced by the motor while the shaft is stationary, the motor is saidto be stalling. This value is limited by the maximum amount of currentthat the amplifier can supply to the motor. What is considered a fastrotation in most industrial applications is at least an order ofmagnitude larger than what one would consider a rapid angular jointvelocity. In industrial applications, motors may move on the order of10³ RPM. In the present invention, the wrist joint is not expected to bemoved to speeds above 25 rads/s or 238 RPM. As shown on the torque-speedcurve, the motor can produce close to the maximal stall torque throughthe entire range of speeds that the present invention is operating atand, thus, the torque-speed relation is not of concern. Rather theconcern is that the maximum toque of 4.8 Nm is sufficient for thepurpose of the present invention. To determine this, the magnitudes oftorques that the motor has to oppose must be analyzed. By Newton'ssecond law,${\sum\tau} = {{\tau_{m} + \tau_{1} + \tau_{f}} = {J_{t}\overset{¨}{\theta}}}$

Rearranging to:$\tau_{m} = {{J_{t}\overset{¨}{\theta}} - \tau_{1} - \tau_{f}}$

wherein τ_(m) is the torque generated by the motor, τ_(l) is the loadtorque applied by the passive and active viscoelastic properties of themuscles and τ_(f) is the negligible friction torque of the mechanicalapperatus. J_(t) is the total inertia of the mechanical apperatus (armrest) J_(m) plus the inertia of the human hand, J_(l). The signs ofτ_(l) and τ_(f) are opposite that of τ_(m) since they oppose the motortorque and rotation direction, thus:$\tau_{m} = {{J_{t}\overset{¨}{\theta}} - \tau_{1} - \tau_{f}}$

The majority of the torque that the motor is required to produce is toovercome the total inertia of the system, J_(t). The inertia of thetypical hand was crudely estimated to be around 0.0015 kg m². Theinertia of the armrest is around 0.005 kg m². J_(t)=J_(m)+J_(l)≈0.0065kg m². Since it is not expected that the wrist joint will be perturbedat accelerations that exceed ±500 rads/s², the maximum torque needed toovercome inertia is${J_{t}{\overset{¨}{\theta}}_{\max}} = {{\left( {0.0065\quad {kgm}^{2}} \right)\left( {500\quad {rads}\text{/}s^{2}} \right)} = {3.25\quad {Nm}}}$

With the load torque applied in response to the perturbations, the motorselected is sufficient for the present invention because the maximumpossible torque that the motor can produce is not expected to beexceeded. Further, based on these specifications, other suitable motorsmay be readily determined by one of skill in the art. The torquecapibility of the motor is expressed in terms of the peak stall torqueand the continuous stall torque. The peak stall torque is the absolutemaximum torque that the motor can produce. This peak torque can only bemaintained for a few seconds before the torque settles down to to thecontinuous stall torque.

The back emf, V_(emf) is an induced voltage by the rotating armaturethat opposes the applied voltage that is determined by: $\begin{matrix}{V_{emf} = {{K_{e}\overset{.}{\theta}\quad {or}\quad {V_{emf}(s)}} = {{K_{e}(s)}\overset{.}{\quad \Theta}\quad \text{in the Laplace domain.}}}} & \text{[eq.~~8]}\end{matrix}$

K_(e) is the emf constant of the motor and θ is the angular velocity ofthe rotor. If the power supplied to the motor is a voltage source andthe voltage supplied at a particular time is V_(m), then:$\begin{matrix}{{V_{m} = {{{L\quad {\overset{.}{I}}_{m}} + {RI}_{m} + V_{emf}} = {{L\quad {\overset{.}{I}}_{m}} + {RI}_{m} + {K_{e}\overset{.}{\theta}\quad {or}}}}}{{V_{m}(s)} = {{\left( {{Ls} + R} \right){I_{m}(s)}} + {{V_{emf}(s)}\quad {in}\quad {the}\quad {Laplace}\quad {domain}}}}} & \left\lbrack {{eq}.\quad 9} \right\rbrack\end{matrix}$

I_(m) is the resulting armature current that can be derived in terms ofV_(m), L and R. L is the inductance, and R is the motor resistance. I isthe rate of change of the armature current. The torque produced by themotor is proportional to the current supplied to it:

τ_(m) =K _(t) I _(m)

or

T _(m)(s)=K _(t) I _(m)(s)

in the Laplace domain  [eq. 10]

where K_(t) is the torque constant of the motor.

Equation 10 can be rewritten as: $\begin{matrix}{{\tau_{m} + \tau_{1}} = {{{B_{m}\overset{.}{\theta}} + {J_{t}\overset{¨}{\theta}\quad \text{or}{T_{m}(s)}} + {T_{1}(s)}} = {\left( {{B_{m}s} + {J_{t}s^{2}}} \right){\Theta (s)}\text{in the Laplace domain.}}}} & \text{[eq.~~11]}\end{matrix}$

The amplifier is the variable power supply to the motor. It amplifiesthe command voltage from the controller into either a voltage or acurrent to the motor 40. The input to the amplifier is the commandvoltage given by the Galil motion controller. There are three modes ofoperation of an amplifier: voltage source mode, current source mode, andvelocity loop mode, all of which are described below:

Voltage Source Mode: In this mode, the output to the motor is a voltage,V_(m), that is some unitless gain, K_(V), times the command voltage,V_(cmd).

V _(m) =K _(V) V _(cmd)  [eq. 12]

The resulting armature current can be derived from equation 9. Theresult is quite cumbersome since equation 9 is a first orderdifferential equation in terms of I_(m) and also depends on the timevarying velocity and the back emf constant, as well as the impedance andresistance of the motor. The resulting motor toque, τ_(m) is simplyproportional to the current I_(m) by the torque constant. It is not easyto derive the torque output of the motor given the amplifier output ofV_(m) because of the complexity of equation 9. The transfer functionthat is the ratio of motor position to voltage command V_(cmd) is:$\begin{matrix}{\frac{\Theta (s)}{V_{cmd}(s)} = \frac{K_{v}K_{t}}{{J_{t}{Ls}^{3}} + {\left( {{B_{m}L_{m}} + {J_{t}R_{m}}} \right)s^{2}} + \left( {{B_{m}R_{m}} + {K_{t}K_{e}}} \right)}} & \left\lbrack {{eq}.\quad 13} \right\rbrack\end{matrix}$

Current Source Mode: The relationship between the motor voltage, V_(m)and motor torque, τ_(m) is a complicated differential equation withrespect to many variables. In current source mode, advantage can betaken of the simple linear relationship between the motor current andtorque output as described by the torque constant, K_(t) as shown inequation 10. Here the amplifier takes the command voltage, V_(cmd) fromthe Galil controller and amplifies it to the current I_(m) the gainK_(a):

I _(m) =K _(a) V _(cmd)  [eq. 14]

The units of K_(a) are [Amps*RMS/Volt]. Here, the transfer functionbetween the motor position and voltage command is: $\begin{matrix}{\frac{\Theta (s)}{V_{cmd}(s)} = \frac{K_{a}K_{t}}{{J_{t}s^{2}} + {B_{m}s}}} & \left\lbrack {{eq}.\quad 15} \right\rbrack\end{matrix}$

Velocity Loop Mode: The velocity loop mode is the same as the currentmode with an added feedback component. A tachometer that senses velocityconverts the velocity into a voltage by the gain, K_(u), and sends itback to the amplifier. The units of K_(u) are [V*s/count]. In this modethe transfer function between the motor position and command voltage is:$\begin{matrix}{\frac{\Theta (s)}{V_{cmd}(s)} = \frac{K_{a}K_{t}}{{J_{t}s^{2}} + {B_{m}s} + {K_{u}K_{a}K_{t}}}} & \left\lbrack {{eq}.\quad 16} \right\rbrack\end{matrix}$

A preferred amplifier for use in the present invention is the KollmorgenBDS4-1033-0101-202B2, which matches the desired motor characteristics.In an illustrative embodiment, the amplifier is set up in the currentmode. This mode makes the analysis of the control system fairly simple.The preferred power supply used for the amplifier is thePSR4/5-112-0102, which can supply a peak current to the above-describedmotor of about 6 Amps RMS, that lasts for around 2 seconds, and amaximum continuous current of about 3 Amps RMS, as shown by the dashedline in FIG. 12.

After setting up the system, prior to operating the motor 40, the gainof the amplifier and other adjustments are set via tuning pots locatedon the chassis of the amplifier. The command scale pot is used to adjustthe velocity loop gain, K_(u). However, since the velocity loop isdisabled, this is not necessary. Pin 18 on Connector C1 of the amplifieris used to monitor the current output relative to common (pin 17). Theoutput is a voltage proportional to the current output where 1V =0.75Amps RMS. This feature is used to monitor the current supplied by theamplifier. The balance pot is used to make sure that when 0V iscommanded from the Galil controller, the current output of the amplifieris 0 Amps RMS and the motor produces no torque. Equation 14 above isreally of the form:

I _(m) =K _(a)V_(cmd)+α  [eq. 14]

where α is the offset current that is adjustable with the balance pot.With the motor enabled, the input to the amplifier, V_(cmd) is set to 0Vand the balance pot is rotated such that I_(m) read 0 Amps to ensure nooffset current. The current limit pot is turned all the way clockwise tomaximize the maximum amount of current that the amplifier can supply (upto ±6 Amps RMS). The most important adjustment is the stability or gainadjustment. This determines the magnitude of K_(a). This adjustableamplifier gain is set experimentally so that the maximum possible valueof V_(cmd) (10 V) would correspond to the maximum current limit(approximately 6A) of the amplifier to take advantage of the full rangeof currents. Thus, the stability pot is adjusted so that value of K_(a)is set to 0.6 Amps RMS/V.

An experiment was conducted where the voltage output from the Galilcontroller was set from 0V to 10V at 0.5V increments. To measure thetorque output of the motor, τ_(m), the motor shaft was kept stationaryusing the mechanical stops. In this state,$\overset{.}{\theta} = {{0\quad {and}\quad \overset{¨}{\theta}} = 0}$

and

τ_(m)=−τ_(m) from equation 11.

This torque was measured with the torque transducer. The motor torqueand amplifier current was measured at each command voltage to producethe plots shown in FIG. 12.

For each command voltage, the correspontind amplifier currents and motortorques were plotted against each other and K_(t) was determined by alinear fit through the origin as shown in FIG. 13.

The maximum torque the motor can produce is limited by the maximumcurrent supplied by the amplifier as defined by equation 10. Thus thepeak transient torque the motor can produce is 4.86 Nm, which lasts forabout 2 seconds before it drops to about 3 Nm.

The controller is the brain of the feedback system, which consists ofhardware and software. It controls the command voltage, V_(cmd) that issent to the amplifier that controls the output torque of the motor shaftsince τ_(m)=K_(t)K_(a)V_(cmd) and K_(a) and K_(t) are known constants.The controller has two inputs at any given discrete time instant, t_(i):the desired position of the motor shaft, θ_(d)(t_(k)) and the actualposition of the motor shaft, θ(t_(i)). The encoder 24 gives the actualposition of the motor shaft 11 with a resolution accuracy of 0.001degrees, as there are 360,000 counts in one revolution of the motorshaft 11. The controller deals with positions in [counts] since that isthe digital output of the encoder/decoder. The command voltage it sendsto the amplifier is some function of the position error, θ_(e) where:

 θ_(e)(t _(i))=θ_(d)(t _(i))−θ(t _(i))

The desired trajectory, θ(t_(k)) is given as a set of digital pointsdt_(mc) seconds apart. This is the period of one cycle of the feed backloop. A new command voltage is sent to the amplifier every dt_(mc)seconds. Thus the frequency of the desired trajectory fed into thecontroller in Hz is ${fhz}_{m\quad c} = \frac{1}{{dt}_{m\quad c}}$

The points in the desired trajectory are either generated offline oronline. If they are generated offline, the exact points of thetrajectory have to predetermined and fed into the controller as anarray. This is what is done, for example, when the desired trajectoriesare the filtered random trajectories that are generated with MATLAB.Since the desired trajectory is periodic, only one period of thetrajectory needs to be inputted to the controller. If the points in thedesired trajectory are dependent on some output state such as torque orcontinuously change with time, then they have to be generated onlinewithin the loop. When the desired trajectory is a sinusoid, each pointof the desired trajectory is generated online based on the time elapsedfrom the start and the user defined amplitude and frequency. If thepoints are to be generated within the loop, then everything from sensingthe necessary output states to computing the next desired position ortorque has to take place in less that dt_(mc) seconds. This is done sothe sampling frequency of the controller is maintained at 128 Hz,otherwise the frequency of the controller needs to be decreased toincrease the amount of time per loop. At substantially low samplingfrequencies, the trajectory will become “choppy” and will not be smoothas defined earlier.

Preferred components of the controller selected for the presentinvention can be obtained from Galil Motion Control. Several high-levelmotion controllers can be considered. Preferably, the controllers can beinterfaced and, thus, controlled by LabVIEW software, which is thepreferred code used in the present device. Preferably, the device isentirely automated and, thus, intervention of the user is minimal ashe/she only gives certain desired inputs such is file identifier and thedesired trajectories for the device to follow. The components of theGalil controller include the following:

the DMC-1410 single-axis motion controller for the ISA bus.

the ICM-1460 interconnect module.

the WSDK-32 Servo Design Kit for Windows 95,98

the Setup-32 software for Windows 95,98

a 37 pin connector cable from the DMC-1410 to the ICM-1460

The controller uses a high level code, which is attached hereto asAppendix A. Of course, controllers having similar components and thatfunction as required by the present invention could also be used.

FIG. 15 shows one preferred embodiment of a block diagram of componentsof the PD closed loop control system used to follow the desiredtrajectories.

As shown in FIG. 15, a PD Compensator is utilized. A preferred PDCompensator for use with the present invention has the followingspecifications:

Manufacturer: Galil Motion Control

Input: position error, θ_(e) [counts]

Output: digital number, N_(d) from −32767 to 32767

Description: motion control software downloaded to PLC takes desiredtrajectory input and feedback from encoder to control motor output.

Transfer Function:$\frac{N_{d}}{\theta_{e}} = \left( {K_{P} + {K_{D}S}} \right)$

Constants: K_(P)=0.75 (chosen) K_(D)=0.1563 s⁻¹ (chosen)

Further shown in FIG. 15 is a Digital to Analog converter (D/A). Apreferred D/A for use with the present invention has the followingspecifications:

Manufacturer: Galil Motion Control

Input: digital number, N_(d) from −32767 to 32767.

Output: command voltage, V_(cmd) from −10 V to 10 V.

Description: converts the digital number from −32767 to 32767 to ananalog voltage output from −10 V to 10 V.

Transfer Function: $\frac{V_{cmd}}{\theta_{e}} = G_{{dac}\quad}$

Constants: G_(dac)=20/65535 (fixed, set by PWM & hardware)

Further shown in FIG. 15 is an amplifier. A preferred amplifier for usewith the present invention has the following specifications:

Manufacturer: Kollmorgen Industrial Drives

Input: command voltage, V_(cmd) from −10 V to 10 V.

Output: motor armature current, i_(m), from −6.16 Amps RMS to 6.16 AmpsRMS.

Description: produces a current proportional to the input voltage. Theadjustable amplifier gain was set experimentally so that the maximumpossible value of V_(cmd) would correspond to the maximum current limitof the amplifier to take advantage of the full range of currents.

Transfer Function: $\frac{I_{m}}{V_{cmd}} = K_{a}$

Constants: K_(a)=0.603 Amps RMS/V (experimentally set)

Further shown in FIG. 15 is a servo motor. A preferred servo motor foruse with the present invention has the following specifications:

Manufacturer: Kollmorgen Industrial Drives

Input: motor armature current, i_(m), from −6.16 Amps RMS to 6.16 AmpsRMS.

Output: torque produced by motor from −5.19 Nm to 5.19 Nm.

Description: produces a torque proportional to the input current.

Transfer Function: $\frac{T_{m}}{I_{m}} = K_{t}$

Constants: K_(t)=0.843 Nm/Amps RMS (fixed, experimentally determined)

Further shown in FIG. 15 is a plant. A preferred plant for use with thepresent invention has the following specifications:

Input: motor torque, T_(m) plus external load torque, T_(l) applied byhuman [Nm].

Output: actual motor shaft position, θ [counts].

Description: the resulting motor shaft position is determined by theinput torques, combined viscosities and inertia of the motor shaft, allmechanical loads attached to the motor shaft including the arm restB_(m) and J_(m). The position is also affected by the inertia of thehand, J_(l) if attached.

Transfer Function:$\frac{\theta}{\left( {T_{l} + T_{m}} \right)} = {\frac{1}{{B_{m}s} + {\left( {J_{l} + J_{m}} \right)s^{2}}} \cdot \frac{4N_{cypr}}{2\pi}}$

Note: all units in the control system are expressed in SI units exceptfor the angular positions which are not in radians but in counts.(4N_(cypr) counts=2π radians)

Constants:

B_(m)=0.003 Nm s/rad experimentally determined, J_(m)≈0.0048 kg m². Theexact value slightly varies depending on the position of the armrest onthe sliding bar.

J_(l)≈0.0023 kg m². The exact value varies depending on the massdistribution, and geometry of the hand of the subject being tested.

N_(cypr)=90,000. This is the number of square wave cycles per revolutionof the incremental encoder. With quadrature edge detection, there are4N_(cypr) counts per revolution.

The entire closed loop PD system preferably has the followingspecifications:

Input: desired motor shaft positions (trajectory), θ_(d) [counts].

Output: actual motor shaft position, θ [counts].

Description: the system implemented is a Linear Time Invariant secondorder model of the system.

Transfer Function: $\begin{matrix}{\frac{\theta}{\theta_{d}}\quad = \quad \frac{{\left( {K_{D}\quad G_{dac}\quad K_{a}\quad K_{t}} \right)\quad s}\quad + \quad {K_{P}\quad G_{dac}\quad K_{a}\quad K_{t}}}{{\frac{\left( {J_{l}\quad + \quad J_{m}} \right)\quad \pi}{2\quad N_{cypr}}\quad s^{2}}\quad + \quad {\left( {{K_{D}\quad G_{dac}\quad K_{a}\quad K_{t}}\quad + \quad \frac{B_{m}\quad \pi}{2\quad N_{cypr}}} \right)\quad s}\quad + \quad {K_{P}\quad G_{dac}\quad K_{a}\quad K_{t}}}} \\{{= \quad}\frac{{24.256\quad s}\quad + \quad 465.705}{{{0.1239\quad s^{2}}\quad + \quad {24.308\quad s}\quad + \quad 465.705}\quad}}\end{matrix}$

It is important that the system characteristic equation or denominatorof R(s), gives poles at: $\frac{\begin{matrix}{{- \left( {{K_{D}G_{dac}K_{a}K_{t}} + \frac{B_{m}\pi}{2N_{cypr}}} \right)} \pm} \\\sqrt{\left( {{K_{D}G_{dac}K_{a}K_{t}} + \frac{B_{m}\pi}{2N_{cypr}}} \right)^{2} - \frac{2{\pi \left( {J_{l} + J_{m}} \right)}\left( {K_{P}G_{dac}K_{a}K_{t}} \right)}{N_{cypr}}}\end{matrix}}{\frac{\left( {J_{l} + J_{m}} \right)\pi}{N_{cypr}}} - {87.56\quad {and}}\quad - 10.73$

after substitution and simplification. Since the roots have negativereal parts, it can be said that the linear time-invariant system isbounded-input, bounded-output and asymptotically stable. The fact thatroots are unequal and have no imaginary parts says that the system isoverdamped (ζ>1). In fact, for the inertia constants used above, ζ=1.6.For smaller inertias, ζ increases slightly and for larger hand inertias,ζ decreases slightly.

Using typical values for constants J_(m), and J_(l) shown above andsubstituting jω for s, the analytical magnitude and phase of thefrequency response is shown in FIG. 16. The analytical step response isalso shown in FIG. 16.

The tuning parameters, K_(p) and K_(d), are chosen so that the systemwill have a fast response time without making them too large. The torquesaturates when small position and velocity error gives rise to a largedigital output N_(d) and a command voltage, V_(cmd) close to +/−10V.Preferably, from 0 to 8 Hz, the magnitude of frequency remains closeto 1. The actual trajectory of the shaft was checked for all types ofdesired trajectories both analytically using the MATLAB function mypid.m(see Appendix A) as well as experimentally by measuring the actualoutput position of the shaft.

A preferred test setup and procedure in accordance with the presentinvention is as follows. The patient being tested with the device of thepresent invention is preferably seated in front of the device. Twosurface electrodes are placed on the patient's arm to monitor activityof the flexors and extensors. The ground electrode is placed on aninactive or neutral site. After the absolute position of the motor shaftis calibrated, the seat is adjusted such that the patent's left or rightarm can be placed in the arm rest 5. At this point, the patient's elbowis preferably at about 145 to 160 degrees. The patient's hand is placedon the hand rest 4 in the neutral supination/pronation position suchthat his or her knuckles are slightly passed between the padded verticalbars 6. The vertical bars 7 are moved closer together and tightened suchthat the hand is comfortably secured in place. The hand rest 4 is thenslid forward or backward so that the wrist joint is lined up with themotor shaft 11. Now the subject is simply told to relax and not toresist the robot's movement as the robot moves to the initial biasposition and follows the desired trajectories as described above. Therobot automatically cycles through the desired trajectories (trials1-10) while collecting torque, position and EMG data. If an additionaltest is to be performed in a different bias position, the patient's armis preferably removed from the apparatus for 2 or 3 minutes before theprocedure is repeated at the new initial bias position.

The bias position is the absolute initial angular position of the wristjoint that the motor shaft 11 moves to before each perturbationtrajectory or trial is carried out. A bias position of 0 degreescorresponds to the neutral flexion/extension position. It is also themedian angular position of the trajectory. The user selects the bias orinitial position, and the motor shaft 11 automatically moves to thatposition before the desired trajectory is carried out. In testingcontrol subjects, selecting a 0-degree bias position is not a problem ascontrol subjects have a wide range of motion. Almost all patients withspasticity and others with neurological problems, however, have alimited range of motion. Their resting position is in a hyper-flexedposition and it is painful or impossible for them to extend their handbeyond the 0-degree absolute position due to the contracted state of theflexor muscles. Due to this constraint, subjects in the patientpopulation were tested in the 30 degree flexed position. Controlsubjects were tested at both 0 and 30 degree flexed bias position. Sincethe left or right arm may be tested, the four possible combinations ofbias positions are shown in FIG. 17.

For more detailed instructions about the procedure, please refer to theuser's manual attached hereto as Appendix B.

Next, the raw data is taken and transformed into the torque offset,elastic stiffness, viscosity and inertia parameters described above andrepresented by the variables τ_(off), K_(H), C_(H), and I_(T)respectively. τ_(off) is the offset torque, K_(H) is the combinedangular elastic stiffness of all the muscles acting on the wrist joint,C_(H) is the combined angular viscosity of all the muscles acting on thewrist joint, and I_(T) is the moment of inertia of the moving hand aswell as the moving mechanical parts of the apparatus (hand rest).

Data obtained from subject JF (a 28 year old, male control subject) isidentified as experiment DMJF. In this experiment, the 10 trial randomtrajectories were run on the right hand at a bias of 0 degrees. Theentire data analysis procedure will be described from obtaining the rawdata until the computation of impendence parameters in equation 1.

As mentioned above, for each perturbation trajectory given, fivemeasurements are collected and stored using LabVIEW and a Nationalinstrument PCI-E series data acquisition board at 500 Hz. A reading fromeach sensor is synchronously collected every 2 ms. These raw data pointsare stored as an n×5 array in ASCII format. Data is collected for 20seconds (n=10,000) per trial. The first data column is just a timestampof when data from each sensor is collected (i.e. 0 ms, 2 ms, 4 ms . . .). The second data column is a set of integer counts from the encoder24, which make up the angular position signal of the motor shaft 11 andthe wrist joint. The third data column is a list of torque readings fromthe torque transducer in Nm. The fourth and fifth data columns areflexor and extensor EMG readings respectively, in mV. Since eachexperiment consists of ten filtered random trajectories, there are tenarrays of raw data stored as DMJFI1t1, DMJFI1t2, . . . , DMJFI1t10.

The raw signals of each of the 10 trials in experiment DMJF are shown inFIGS. 18-27.

In one preferred embodiment, the device utilizes a filter. One preferredfilter is a Savitsky-Golay filter (Savitsky & Golay, Press et. al:Numerical Recipes). A filter can provide a number of advantages. Thediscretization of position gives a finite resolution of position basedon the resolution of the encoder. In the present invention, thisresolution is 0.001 degrees. Filtering enables displacement positions tofall between these increments. Since the resolution is very high, thisis particularly important when using low-resolution encoders, with a lowncypr. Further, filters can assist in the timing of data acquisition.Still further, filters are useful in obtaining smooth first and higherorder derivatives of a discrete signal. A crude method to obtain thederivative (slope or tangent) at a point in a discrete signal is to takethe difference between neighboring points and then divide it by thesampling period. This is based on the assumption that the twoneighboring points are sufficiently close enough such that the points inbetween lie on a straight-line between the two. Sharp discontinuities inthe derivative signal are observed using this technique. As in mostsmooth signals, the present displacement signal between any two pointsis actually a curve. Thus, the filter utilizes a better technique toobtain derivatives by fitting a polynomial through a window ofneighboring points. The Savitsky-Golay filter, for example, is asmoothing filter that is useful when determining derivatives such asvelocity and acceleration offline when the raw displacement signal hasalready been collected and stored. This smoothing technique fits apolynomial to a moving window of data using least squares fitting.Derivatives, as well as the function value, can be determined at thecenter point of each window. Since the coefficients of the fittedpolynomial are linear in the data values, the filter coefficients can bepre-computed. The filter coefficients are then convoluted with the rawdigital signal to obtain the actual smoothed version of the raw waveformas well as well as its derivatives at each center point.

The Savitsky-Golay algorithm operates as follows with the followingvariables:

d_(raw): the column vector of actual displacement data points

N: total number of points in the raw digital displacement signal (rads)

ord: the order of the polynomial.

wins: the number of equally spaced points to be fit by the polynomial,the size of the window. This number is always odd.

y_(i): a sub vector of d_(raw) of length wins with i being the index ofthe center point.

hw=(wins−1)/2: half width of the window.

The matrix X of dimension wins by (ord+1) is taken as$X_{{winsXord} + 1} = {\begin{matrix}{- {hw}^{0}} & {- {hw}^{1}} & \cdots & {- {hw}^{{ord} - 1}} & {- {hw}^{ord}} \\{- \left( {{hw} - 1} \right)^{0}} & {- \left( {{hw} - 1} \right)^{1}} & \cdots & {- \left( {{hw} - 1} \right)^{{ord} - 1}} & {- \left( {{hw} - 1} \right)^{ord}} \\\vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & 0 \\\left( {{hw} - 1} \right)^{0} & \left( {{hw} - 1} \right)^{1} & \cdots & \left( {{hw} - 1} \right)^{{ord} - 1} & \left( {{hw} - 1} \right)^{ord} \\{- {hw}^{0}} & {- {hw}^{1}} & \cdots & {- {hw}^{{ord} - 1}} & {- {hw}^{ord}}\end{matrix}}$

The sum of the squared error, S, is minimized between the elements ofy_(i) and y_(i), an estimate of y_(i):

y _(i) =X*F*y _(i)

where F is a (ord+1) by wins matrix. The sum of the squared errors canbe written as:

S=(y _(i) −y _(i))^(T)*(y _(i) −y _(i))=(y _(i) −X*F*y _(i))^(T)(y _(i)−X*F*y _(i))

A matrix F is evaluated such that the least possible sum-squared error,S, is obtained. To do this, the derivative of S is set with respect to Fto 0, dS/dF=0. This yields the pseudo-inverse solution,

F=(X ^(T) *X)⁻¹ *X ^(T)

The columns of F are then reversed and the k^(th) derivative is thengiven by k! multiplied with the (k+1)^(th) row of F divided by(dt_daq)^(k) convolved with the raw displacement data. Before applyingthe Savitsky-Golay filter, the encoder count signal is converted intoradians. Since the resolution of the encoder gives us 1000 counts for adegree of angle, the raw encoder count signal is multiplied by$\frac{\pi}{180000}$

to give the raw displacement data in radians in the column vector,rawrads. In equation form, the i^(th) element of the k^(th) derivativeusing the Savitsky-Golay filter is given as:${\frac{^{k}\left( {r\quad a\quad w\quad r\quad a\quad d\quad s} \right)}{t^{k}}(i)} = {\left( \frac{1}{{{t\_}}{a}\quad q^{k}} \right)*{k!}*{\sum\limits_{j = {{- h}\quad w}}^{h\quad w}\left( {r\quad a\quad w\quad r\quad a\quad d\quad s_{i}*F_{{k + 1},{i + 1 - j}}} \right)}}$

A third order polynomial (ords=3) is used with a window size of 21points. Thus, F is a 4 by 21 matrix. From the equation above, the i^(th)element of the filtered displacement signal is given as:${d\quad i\quad s\quad {p(i)}} = {\sum\limits_{j = {- 10}}^{10}\left( {r\quad a\quad w\quad r\quad a\quad d\quad s_{i}*F_{1,{i + 1 - j}}} \right)}$

the velocity signal is given as:${{vel}(i)} = {{\frac{\left( {r\quad a\quad w\quad r\quad a\quad d\quad s} \right)}{t}(i)} = {\left( \frac{1}{0.002} \right)*{\sum\limits_{j = {- 10}}^{10}\left( {r\quad a\quad w\quad r\quad a\quad d\quad s_{i}*F_{2,{i + 1 - j}}} \right)}}}$

and the acceleration signal is given as:${{acc}(i)} = {{\frac{^{2}\left( {r\quad a\quad w\quad r\quad a\quad d\quad s} \right)}{t^{2}}(i)} = {\left( \frac{1}{4*10^{- 6}} \right)*2*{\sum\limits_{j = {- 10}}^{10}{\left( {r\quad a\quad w\quad r\quad a\quad d\quad s_{i}*F_{3,{i + 1 - j}}} \right).}}}}$

Since the units of the elements in the raw rads vector are in rads orradians, elements in the displacement, velocity, and acceleration arrayare in rads, rads/s, and rads/s² respectively. The beginning and the endof each vector are truncated by a number of points equal to the halfwidth of the window. The frequency response of the Savitsky-Golay filterusing a third order polynomial with a window size of 21 points (0.04 s)in the frequency domain can be seen using the freqz MATLAB functionwhich shows the Z-transform digital filter frequency response asfollows:

>>sav=sav_golay(3,21);

>>[h,f]=freqz(sav(1,:),1);

>>plot(f*500/(2*pi),abs(h),‘k’);

This produces the frequency response plot shown in FIG. 28. There is nophase lag since the window is symmetric with respect to the point atwhich the polynomial is evaluated (center point).

Unlike digital signals, such as that from the encoder 24, analog signalsfrom torque transducers have high frequency noise components that do notgive a true representation of the actual torque exerted. Low passdigital filters are often used offline to attenuate the high frequenciesin the torque signal. A zero lag 8^(th) order digital Butterworth filtercan be used to greatly attenuate these high frequencies.

In MATLAB, [b,a]=butter(n,Wn) designs an order n lowpass digitalButterworth filter with cutoff frequency Wn. It returns the filtercoefficients of length n+1 in row vectors b and a, with coefficients indescending powers of z:${H(z)} = {\frac{B(z)}{A(z)} = \frac{{b(1)} + {{b(2)}z^{- 1}} + \ldots + {{b\left( {n + 1} \right)}z^{- n}}}{1 + {{a(2)}z^{- 1}} + \ldots + {{a\left( {n + 1} \right)}z^{- n}}}}$

The cutoff frequency is that frequency where the magnitude response ofthe filter is {square root over ({fraction (1/2)})}. For example, forMATLAB's butter function, the cutoff frequency, Wn must be a numberbetween 0 and 1, where 1 corresponds to half the sampling frequency (theNyquist frequency). A preferred embodiment of the present invention usesa Wn of 0.075, which corresponds to a cutoff frequency of 18.75 Hz (seegetstatesgbw.m in Appendix A). Normally, Butterworth filters have aphase response where there is a lag between the filtered and the rawsignal as a function of frequency of the raw signal. In other words thefiltered torque signal is nonlinearly shifted in time with respect tothe original signal as a function of frequency. Thus, the filteredtorque signal also lags behind the encoder signal. These phase lags areparticularly significant in high order filters. This lag is problematicbecause it is important for the torque waveform to be synchronized withthe position velocity and acceleration waveforms in order to evaluatethe parameters in the viscoelastic model described herein. However, thefiltfilt function in MATLAB is prefereably used. y=filtfilt(b,a,x)performs zero-phase digital filtering by processing the input data inboth the forward and reverse directions.

After filtering in the forward direction, the filtered sequence isreversed and it runs it back through the filter. The resulting sequencehas precisely zero-phase distortion and double the filter order (A. V.Oppenheim, and R. W. Schafer. Discrete-Time Signal Processing. EnglewoodCliffs, N.J.: Prentice Hall, 1989. Pgs. 311-312).

Stiffness, viscosity, and inertial parameters are obtained from positionand torque signals obtained from the optical encoder and torquetransducer. Velocity and acceleration signals are also required, asshown in equation 1. These signals are derived inherently from theposition signal via the Savitsky-Golay filter described above. Filteringof the flexor and extensor EMG signals is not crucial since thesesignals are not necessary for post processing. However, monitoring liveEMG information can be used to monitor possible voluntary contractions.Ideally the muscle should be fully relaxed before perturbing the wrist,since spasticity is being measured, which is an involuntary reflex.Thus, the EMG information can be used live to ensure proper integrity ofthe data collected. Since flexor and extensor EMG activity is monitored,the EMG data can also be collected and stored. Post processing on EMGsignals can be useful to answer more basic questions, such as theinvestigating reflex delay times from the onset of single or repetitivemuscle perturbations to muscle contraction. Thus, the total time ittakes for the entire reflex loop can be measured. Researchers have shownthat two separate response times can be measured, one due to spinalpathways to the spinal cord and back, and one due to supra-spinalpathways which involve upper motor neurons from the brain (See J. C.Houk: Participation of Reflex Mechanisms and Reaction-Time Processes inthe Compensatory Adjustments to Mechanical Disturbances. Cerebral MotorControl in Man: Long Loop mechanisms. Prog. clin. Neurophysiol., vol 4,Ed. J. E. Desmedt, pp 193-215, Krager Basel, 1978; W. G. Tafton, P.Bawa, I. C. Bruce, R. G. Lee: Long Loop Reflexes in Monkeys: AnInterpretative Base for Human Reflexes. Cerebral Motor Control in Man:Long Loop mechanisms. Prog. clin. Neurophysiol., vol 4, Ed. J. E.Desmedt, pp 229-245, Krager Basel, 1978). Naturally, the supra-spinalresponse occurs after the spinal response.

EMG signals are extremely noisy, particularly if surface electrodes areused. EMG signals measured with surface electrodes are influenced byhair, oil, lotions and dead skin. The actual action potentials sent tothe muscles are internally amplified by 1000, typically resulting in asignal from 0 to ±5 V. Preferably, internal to the Delsys 2-channel EMGunit is a 20-450 Hz band pass filter. Once the data is collected throughthe data acquisition board and stored, further filtering is thenperformed offline in MATLAB as follows.

First a moving root-mean-square (RMS) window is applied to the rawoutput of the EMG signal where the size of the window is 15 points or 28milliseconds. Thus, the k^(th) sample point is of the RMS EMG signal isequal to${E\quad M\quad {G_{R\quad M\quad S}(k)}} = \sqrt{\frac{1}{15}{\sum\limits_{i = {k - 7}}^{k + 7}\left( {E\quad M\quad {G_{R\quad A\quad W}(k)}} \right)^{2}}}$

After the RMS EMG signal is obtained, it is passed through a low passbutterworth filter as applied to the torque signal. The result gives thefinal filtered EMG signal. The cutoff frequency used in this butterworthfilter is 37.5 Hz.

The displacement, velocity and acceleration signals using theSavitsky-Golay filter on the encoder signal for each of the 10 trialsare shown in FIGS. 30-29. The filtered torque and EMG signals are alsoshown. Filtered data obtained from two cycles or periods of each randomdisplacement signal are used for further processing.

Now all the proper state measurements$\left( {{\tau_{s}\left( t_{i} \right)},{\theta \left( t_{i} \right)},{\overset{.}{\theta}\left( t_{i} \right)},{{and}\quad {\overset{¨}{\theta}\left( t_{i} \right)}}} \right)$

are known after the signal processing described above, and the data cannow be fit to the model of the present invention:${\tau_{s}\left( t_{i} \right)} = {\tau_{off} + {K_{H}*{\theta \left( t_{i} \right)}} + {B_{H}*{\overset{.}{\theta}\left( t_{i} \right)}} + {J_{T}*{\overset{¨}{\theta}\left( t_{i} \right)}}}$

where τ_(S) is the filtered torque signal from the transducer. K_(H) andB_(H) are the angular stiffness and viscosity of the combined flexor andextensor muscle groups that act on the wrist joint. J_(T) is thecombined inertia of oscillating appendage, J_(H), and the rotatingcomponents of the apparatus, J_(A). J_(T)=J_(H)+J_(A). τ_(off), theoffset torque is equal to −K_(H)*θ_(RP) as described above. θ is theangular displacement of the system.$\overset{.}{\theta}\quad {and}\quad \overset{¨}{\theta}$

are velocity and acceleration, respectively which are the first andsecond derivatives of displacement. The data was fit using a linearleast squares fit pseudo-inverse (See A. Bjorck, Numerical Methods forLeast Squares Problems, SIAM, Philadelphia, 1996; H. Anton, C. Rorres,Elementary Linear Algebra, Applications version, sixth edition, Wiley &Sons, New York, 1991), as follows. For a set of discrete samples, thefollowing matrices are constructed:${T\quad s} = {\left| \begin{matrix}{\tau \quad {s\left( t_{1} \right)}} \\{\tau \quad {s\left( t_{2} \right)}} \\{\tau \quad {s\left( t_{3} \right)}} \\\vdots \\{\tau \quad {s\left( t_{n} \right)}}\end{matrix} \middle| {}_{n \times 1}\quad {{and}\quad \Psi} \right. = \left| \begin{matrix}1 & {\theta \left( t_{1} \right)} & {\overset{.}{\theta}\left( t_{1} \right)} & {\overset{¨}{\theta}\left( t_{1} \right)} \\1 & {\theta \left( t_{2} \right)} & {\overset{.}{\theta}\left( t_{2} \right)} & {\overset{¨}{\theta}\left( t_{2} \right)} \\1 & {\theta \left( t_{3} \right)} & {\overset{.}{\theta}\left( t_{3} \right)} & {\overset{¨}{\theta}\left( t_{3} \right)} \\\vdots & \vdots & \vdots & \vdots \\1 & {\theta \left( t_{n} \right)} & {\overset{.}{\theta}\left( t_{n} \right)} & {\overset{¨}{\theta}\left( t_{n} \right)}\end{matrix} \right|_{n \times 4}}$

It is possible to obtain the set of three parameters,$\overset{\sim}{N} = {\begin{matrix}\begin{matrix}\begin{matrix}{\overset{\overset{\sim}{\hat{}}}{o}}_{off} \\{\overset{\sim}{K}}_{H}\end{matrix} \\{\overset{\sim}{B}}_{H}\end{matrix} \\{\overset{\sim}{J}}_{T}\end{matrix}}_{4 \times 1}$

that best fits the present model in equation 1:${\tau_{S}\left( t_{i} \right)} = {\tau_{off} + {K_{H}*{\theta \left( t_{i} \right)}} + {B_{H}*{\overset{.}{\theta}\left( t_{i} \right)}} + {J_{T}*{\overset{¨}{\theta}\left( t_{i} \right)}}}$

using the pseudo-inverse which minimizes the sum squared error, apositive scalar value. In matrix equation form the sum-squared error canbe written as:

SSE=(Ô _(S) −Ø*Ñ)^(T)*(Ô _(S) —Ø*Ñ)

The matrix P that minimizes SSE is evaluated as:

Ñ=pinv(Ø)*Ô _(S)=(Ø^(T)*Ø)⁻¹*Ø^(T) *Ô _(S)  [eq. 15]

First the pseudo-inverse is performed to obtain one set of parametersutilizing all data over 2 whole periods (15.625 s or n=7815) thatconstruct one large Ô_(S) (7815×1) vector of torque measurements and onelarge Ø (7815×4) state matrix. Let us denote the four evaluatedparameters as {tilde over (τ)}_(off), {tilde over (K)}_(H), {tilde over(B)}_(H), and {tilde over (J)}_(T). These evaluated parameters for theexperiment DMJF whose data is shown below.

TABLE 2 The use of overall data in each trial to obtain one set ofbest-fit parameters for each trial for experiment DMJF. {tilde over(τ)}_(off) {tilde over (K)}_(H) {tilde over (B)}_(H) {tilde over(J)}_(T) Torque Offset Stiffness Viscosity Total Inertia Trial [N * m][N * m/rad] [N * m * s/rad] [kg * m²]  1 −0.03667 0.76224 0.050090.00687  2 −0.03744 0.76663 0.04372 0.00673  3 −0.03740 0.81638 0.045030.00682  4 −0.04743 0.89976 0.04170 0.00701  5 −0.04922 0.83478 0.043550.00701  6 −0.05149 0.88837 0.04229 0.00707  7 −0.06263 0.92421 0.045940.00720  8 −0.06939 0.89970 0.04502 0.00715  9 −0.06859 0.94812 0.050220.00730 10 −0.07785 1.01639 0.05268 0.00737 Average −0.05381 0.875660.04603 0.00705 Standard   0.01496 0.08074 0.00372 0.00021 dev.

To determine whether the model used in the present invention to obtainoverall parameters for each trial is a reasonable one, the residualtorque, τ_(R) is plotted. The residual torque is the amount of torquethat is unaccounted for by the model and for a particular time sample.In equation form, it is defined as:${\tau_{R}\left( t_{i} \right)} = {{\tau_{s}\left( t_{i} \right)} - \left\lbrack {{\overset{\sim}{\tau}}_{off} + {{\overset{\sim}{K}}_{H}*{\theta \left( t_{i} \right)}} + {{\overset{\sim}{B}}_{H}*{\overset{.}{\theta}\left( t_{i} \right)}} + {{\overset{\sim}{J}}_{T}*\overset{¨}{\theta}\left( t_{i} \right)}} \right\rbrack}$

From the residual torque plots in FIG. 40, it is shown that the amountof measured torque, õ_(S)(t_(i)) that is not accounted by the model (theresidual, τ_(R)(t)) is only about 4% on average. This demonstrates thatthe present model described by equation 1 fits to the data well, evenwhen estimating only one set of parameters, {tilde over (τ)}_(off),{tilde over (K)}_(H), {tilde over (B)}_(H) and {tilde over (J)}_(T) overthe entire data set in each trial. If, however the residuals are large,this would indicate that the equation of the present model isinadequate, or that a single set of time invariant, four parametersapplied to the state data$\left( {{\theta (t)},{\overset{.}{\theta}(t)},{{and}\quad {\overset{¨}{\theta}(t)}}} \right)$

over the whole trial is insufficient to describe the torque output. Inother words, if the residual is large using only one set of timeinvariant parameters, the model may be good, but the impedanceparameters of torque offset, stiffness and viscosity may be changing intime, within a given trial. To estimate how the impedance parameterschange within a given trial, rather than have a single T_(S), P_(S) andΨ matrix over the whole trial as shown above, multiple T_(S), P_(S) andΨ matrices can be used by applying a smaller, one second moving window(501 samples) over the data as follows. First it is assumed that thetotal inertia does not change within each trial. This is a reasonableassumption since the distribution of mass of the moving parts about thecenter of the motor shaft does not change and the subject's hand issecured in place with the vertical bars 6. For each trial, {tilde over(J)}_(T) is used, which is evaluated above as the total inertia of thesystem. It is possible for the other three parameters to change within atrial. For a particular point in time, these now time varying parametersare denoted with subscripts as τ_(off)(t_(i)), K_(H)(t_(i)),B_(H)(t_(i)). Using these time varying parameters and subtracting thetorque due to inertia results in a slightly modified equation of themodel:${{\tau_{S}\left( t_{i} \right)} - {J_{T}*{\overset{¨}{\theta}\left( t_{i} \right)}}} = {{\tau_{off}\left( t_{i} \right)} + {K_{H}*{\theta \left( t_{i} \right)}} + {B_{H}*\overset{.}{\theta}\left( t_{i} \right)}}$

The torque data set with out the contribution of inertia is now fit asfollows: The data is fit using a linear least squares fit pseudo-inversesimilar to above, but this is performed multiple times, once for eachwindow as follows:${{T_{si} - {{\overset{\sim}{J}}_{T}*{\overset{¨}{\Theta}}_{i}}} = {{\begin{matrix}{{\tau_{s}\left( t_{i - 250} \right)} - {{\overset{\sim}{J}}_{T}*{\overset{¨}{\theta}\left( t_{i - 250} \right)}}} \\\vdots \\{{\tau_{s}\left( t_{i} \right)} - {{\overset{\sim}{J}}_{T}*{\overset{¨}{\theta}\left( t_{i} \right)}}} \\\vdots \\{{\tau_{s}\left( t_{i + 250} \right)} - {{\overset{\sim}{J}}_{T}*{\overset{¨}{\theta}\left( t_{i + 250} \right)}}}\end{matrix}}_{501 \times 1}{and}}}\quad$$\Psi_{i} = {\begin{matrix}1 & {\theta \left( t_{i - 250} \right)} & {\overset{.}{\theta}\left( t_{i - 250} \right)} \\\vdots & \vdots & \vdots \\1 & {\theta \left( t_{i} \right)} & {\overset{.}{\theta}\left( t_{i} \right)} \\\vdots & \vdots & \vdots \\1 & {\theta \left( t_{i + 250} \right)} & {\overset{.}{\theta}\left( t_{i + 250} \right)}\end{matrix}}_{501 \times 3}$

where i goes form 251 to 8064 that cover two periods of the randomtrajectory or 15.63 seconds). For each i, the matrices$\left\lbrack {T_{si} - {{\overset{\sim}{J}}_{T}*{\overset{¨}{\Theta}}_{i}}} \right\rbrack$

and Ø_(i) are obtained and P_(i) is computed for each i as follows:${\overset{\sim}{N}}_{i} = {\begin{matrix}{{\hat{\tau}}_{OFF}\left( t_{i} \right)} \\{{\hat{K}}_{H}\left( t_{i} \right)} \\{{\hat{B}}_{H}\left( t_{i} \right)}\end{matrix}}_{3 \times 1}$

where P_(i) is calculated similar to P above as follows:

Ñ _(i) =pinv(Ø_(i))*Ô _(S) _(i) =(Ø_(i) ^(T)*Ø_(i))⁻¹*Ø_(i) ^(T)*Ô_(S)_(i)

Since I goes from 251 to 8064, there are 7813 {circumflex over(τ)}_(off), {circumflex over (K)}_(H), and {circumflex over (B)}_(H)values per trial. The mean and variance or standard deviation of thesevalues for each trial is determined. If the standard deviation of thesevalues: std{{circumflex over (τ)}_(off)(t_(i))}, std{{circumflex over(K)}_(H)(t_(i))}, and std{{circumflex over (B)}_(H)(t_(i))} are small,this indicates that the stiffness and viscosity of the subject's flexorsand extensors remain fairy constant within each trial and the torqueresiduals of the time invariant model would be small as shown in FIG.40. Note also that the mean of the time-variant impedance parametersvalues (denoted with the {circumflex over (x)}) should be very close tothe evaluated time-invariant parameters (denoted with the {tilde over(x)}), that is:

mean{{circumflex over (τ)}_(off)(t_(i))}≈{tilde over (τ)}_(off);

mean{{circumflex over (K)}_(H)(t_(i))}≈{tilde over (K)}_(H); and

mean{{circumflex over (B)}_(H)(t_(i))}≈{tilde over (B)}_(H);

This is shown in FIG. 41.

The first row of bar plots in FIG. 41 show the results for the torqueoffset, stifffiess, and viscosity parameters obtained for each of theten trials in experiment DMJF. The asterisk ‘*’ plots the values of{tilde over (τ)}_(off), {tilde over (K)}_(H), and {tilde over (B)}_(H)where one set of best-fit parameters over the whole trial as describedabove. The location of the tip of each bar graph representsmean{{circumflex over (τ)}_(off)(t_(i))}, mean{{circumflex over(K)}_(H)(t_(i))}, and mean{{circumflex over (B)}_(H)(t_(i))}where a setof {tilde over (τ)}_(off), {circumflex over (K)}_(H), and {circumflexover (B)}_(H) is obtained for each moving window of 501 data samples asdescribed above. The height of the error bar indicates the standarddeviation of these values, std{{circumflex over (τ)}_(off)(t_(i))}std{{circumflex over (K)}_(H)(t_(i))}, and std{{circumflex over(B)}_(H)(t_(i))}.

The eleventh bar shows a summary of all ten trials. The asterisk on theeleventh bar shows the average of the ten {tilde over (τ)}_(off), {tildeover (K)}_(H), or {tilde over (B)}_(H) values, or in equation form:${{\frac{1}{10}{\sum\limits_{k = 1}^{10}\left( {{\overset{\sim}{\tau}}_{off}\quad {for}\quad {trial}\quad k} \right)}},{\frac{1}{10}{\sum\limits_{k = 1}^{10}\left( {{\overset{\sim}{K}}_{H}\quad {for}\quad {trial}\quad k} \right)}}\quad,\quad {and}}\quad$$\quad {\frac{1}{10}{\sum\limits_{k = 1}^{10}\left( {{\overset{\sim}{B}}_{H}\quad {for}\quad {trial}\quad k} \right)}}$

The location of the tip of the eleventh bar is the average of the heightof the ten bar graphs, or in equation form:${{\frac{1}{10}{\sum\limits_{k = 1}^{10}\left( {{mean}\left\{ {\hat{\tau}}_{off} \right\} \quad {for}\quad {trial}\quad k} \right)}},{\frac{1}{10}{\sum\limits_{k = 1}^{10}\left( {{mean}\left\{ {\overset{\sim}{K}}_{H} \right\} \quad {for}\quad {trial}\quad k} \right)}},\quad {{and}\quad \frac{1}{10}{\sum\limits_{k = 1}^{10}\left( {{mean}\left\{ {\hat{B}}_{H} \right\} \quad {for}\quad {trial}\quad k} \right)}}}\quad$

As mentioned earlier, there are 7813 {circumflex over (τ)}_(off),{circumflex over (K)}_(H), and {circumflex over (B)}_(H) values pertrial. If all ten trials are grouped together, then there are 78130{circumflex over (τ)}_(off), {circumflex over (K)}_(H), and {circumflexover (B)}_(H) values. The height of the error bar of the eleventh barplot is the standard deviation of all 78130 values.

Experiment DMJF was performed on the right arm of a control adultsubject tested at a bias of 0 degrees. Each control subject was alsotested at a bias position of degrees flexion. The bottom row of graphsin FIG. 41 shows the results of parameters for two sets of experiments,DMJF, and DNJF, as labeled on the left of the graphs. The ‘M’ in theparentheses, indicates that the subject JF is a male. The bar graphs forexperiment DMJF are plotted in white (bias of 0 degrees) and bar graphsfor experiment DNJF are plotted in grey (bias of 30 degrees flexion) asindicated by the legend.

The results of subject JF show that the stiffness values when tested ata bias of 30 degrees are markedly higher than when tested at a bias of 0degrees in all ten trials. The average {tilde over (K)}_(H) across allten trials are 0.876 N*m/rad and 1.30 N*m/rad at a bias of 0 and 30degrees respectively. The asterisks on the eleventh bars show thesevalues. Thus, the cumulative stiffness of all the muscles acting on thewrist joint is 48% higher at a bias of 30 degrees flexion.

Similarly, the viscosity values are markedly higher when tested at abias of 30 degrees flexion than at a bias of 0 degrees for all tentrials. The average {tilde over (C)}_(H) across all ten trials are0.0460 N*m*s/rad and 0.0714 N*m*s/rad at a bias of 0 and 30 degreesrespectively and the cumulative stiffness of all the muscles acting onthe wrist joint is 55% higher at a bias of 30 degrees flexion.

The resulting torque offsets, {tilde over (τ)}_(off) tell us somethingvery interesting. The first noticeable thing about these values is notonly that they are different in magnitude but also opposite in sign.τ_(off) is an important variable that tells us the angular position thehand prefers to be relative to the selected bias position or origin, θ₀.If the bias position was chosen to be the resting position of the wristjoint, then in the equation below θ_(RP)=0,

τ_(off) =−K _(H)*θ_(RP)

and ideally, τ_(off) would be 0 Nm. For a given angular stiffness, thefurther the bias position is selected from the preferred restingposition of the wrist joint, the greater the torque offset, τ_(off) willbe. The preferred resting position of the wrist joint is determined bythe relative passive stiffness of flexor muscle group to the extensormuscle group. It should also be noted that τ_(OFF) is also dependent onthe stiffness K_(H). Since the active stiffness of muscles areinfluenced when they are perturbed, τ_(OFF) is influenced by these wristperturbations or oscillations. If, however, the bias position is veryclose to the resting position of the wristjoint, τ_(off) will be veryclose to 0 Nm regardless of the magnitude of K_(H) since θ_(RP)≈0. Apositive torque offset indicates that this contribution of torque isapplied in the clockwise direction by the subject's arm and a negativetorque offset indicates that this contribution of toque is applied inthe counter clockwise direction. This information is important wheninterpreting the sign of {tilde over (τ)}off. The subject JF's right armwas tested in experiments DMJF and DNJF. As shown in the plots in FIG.41, when tested at a bias of 0 degrees, the torque offset values arevery slightly negative. The mean overall torque offset, {tilde over(τ)}_(off) for all ten trials shown by the asterisk on the eleventh baris equal to −0.0538 Nm at a bias of 0 degrees, and the contribution oftorque by {tilde over (τ)}_(off) is in the counter clockwise direction.At a bias of 30 degrees flexion, this value is equal to 0.110 Nm and thecontribution of torque by {tilde over (τ)}_(off) in this case is in theclockwise direction. This means that there is less resistance in thewrist flexors at a bias of 0 than of 30 degrees flexion and the restingposition of the hand of the control subject is close to the neutralflexion-extension position. The magnitude of these numbers is small ascompared to the torque offset values in the spastic population as willbe shown. It is important to note that if contra-lateral left arm of thesubject was tested instead, then the sign of τ_(off) is expected toswitch by virtue of symmetry. Thus, in order to compare the left hand tothe right hand, if the left hand of a subject is tested, then the signof τ_(off) is switched. Hence, the sign of τ_(off) no longer representsclockwise or counterclockwise torque. Rather, a positive torque offsetindicates that the contribution of τ_(off) is in the direction towardsthe extension direction of the hand being tested relative to theselected bias position. Conversely, a negative torque offset indicatesthat the contribution of τ_(off) is in the direction towards the flexiondirection of the hand being tested relative to the selected biasposition.

Now, the results of each individual control adult will be shown. Eitherthe left or right hand of 31 control adults was tested, at a biasposition of both 0 and 30 degrees flexion. 10 of the subjects were maleand 21 were female. The average age of control subjects tested was31.29, the youngest was 18 years of age and the eldest was 60 years ofage. The plots of the results for each individual subject are shown inFIGS. 42-52. All graphs that are plots of results of the same variableare plotted on the same scale so that visual comparisons can be made.

As shown in FIGS. 42-52, when looking at the average of {tilde over(τ)}_(off), {tilde over (K)}_(H) and {tilde over (B)}_(H) over all tentrials for each experiment (denoted by the asterisk over the eleventhbar labeled trial 11), the following conclusions can be made. At a biasposition of 0, the magnitude of {tilde over (τ)}off across all tentrials, for all 31 subjects is negligible, some being positive in value,some negative. This indicates that the resting position of a controlsubject's hand is very close to the neutral flexion/extension wheremuscles acting on the wrist exhibit the least amount of tension. At abias of 30 degrees the torque-offset values are significantly higher andpositive in sign. This indicates that the hand exhibits more elastictension at this bias position. The positive sign indicates that the handprefers to be at the neutral 0 degree bias position. All 31 subjects hada higher stiffness at a 30 degree flexion bias than at a 0 degree bias.30 subjects had a higher mean viscosity value at a 30 degree flexionbias than at a 0 degree bias. Only 1 subject (MT) had a lower meanviscosity value at a 30 degree flexion bias than at a 0 degree bias. Theaverage and standard deviation (denoted by the height of the error bars)of {tilde over (τ)}_(off), {tilde over (K)}_(H) and {tilde over (B)}_(H)values across all control adults for each trail are shown in FIG. 53.The eleventh bar shows the overall average and standard deviation acrossall control adults, across all ten trials.

Null-hypothesis significance tests can be used to quantitatively testhow significant these different observations are between testing at thebias of 30 degrees flexion versus the bias of 0 degrees. For eachexperiment, the mean {tilde over (τ)}_(off), {tilde over (K)}_(H), and{tilde over (B)}_(H) are taken over the ten trials and these parametersare separated into two groups according to bias position. Thus, thereare 31 mean{{tilde over (τ)}_(off)}, mean{{tilde over (K)}_(H)}, andmean{{tilde over (B)}_(H)} parameters for each of the two biaspositions. The null hypotheses test indicates that the mean torqueoffsets, mean stiffnesses and mean viscosities are not significantlyhigher at a bias of 30 degrees than they are at a bias of 0 in controlsubjects. From the null hypothesis test, the p-value is obtained,wherein the p-value is the probability of observing the given result bychance given that the null hypothesis is true. Small values of p castdoubt on the validity of the null hypothesis. The probability of thenull hypothesis stating that the average torque offset values are higherat 0 than at 30 degrees is true is p=9.2*10⁻⁴ The probability of thenull hypothesis stating that the average stiffness values are higher at0 than at 30 degrees is true is p=2*10⁻⁵ The probability of the nullhypothesis stating that the average viscosity values are higher at 0than at 30 degrees is true is p=4.2*10⁻⁵. Since a resulting p-value of5*10⁻³ is considered to be very significant, there is no question thatthe increase in all three parameters {tilde over (τ)}_(off), {tilde over(K)}_(H), and {tilde over (B)}_(H) is extremely significant when testingat a bias position of 30 degrees flexion rather than 0 in the controladult subject population.

Table 3 below shows the average and standard deviation of {tilde over(τ)}_(off), {tilde over (K)}_(H) and {tilde over (B)}_(H) values acrossall ten trials across all control adult subjects which are indicated bythe height of the 11^(th) bars shown in FIG. 53. Table 4 shows the rangeof average {tilde over (τ)}_(off), {tilde over (K)}_(H) and {tilde over(B)}_(H) values across all ten trials in the control adult group.

TABLE 3 {tilde over (τ)}_(off) [N * m]: Bias Torque offset {tilde over(K)}_(H) [N * m/rad]: {tilde over (B)}_(H) [N * m * s/rad]: positiontowards extension Stiffness Viscosity  0 degrees −0.0189 ± 0.0609 0.789± 0.309 0.0353 ± 0.0127 30 degrees   0.1667 ± 0.0936 1.305 ± 0.4800.0512 ± 0.0175 flexion

TABLE 4 {tilde over (τ)}_(off) [N * m]: Bias Torque offset {tilde over(K)}_(H) [N * m/rad]: {tilde over (B)}_(H) [N * m * s/rad]: positiontowards extension Stiffness Viscosity  0 degrees −0.123 to 0.089 0.233to 1.347 0.021 to 0.059 30 degrees 0.0241 to 0.361 0.482 to 2.345 0.0274to 0.105  flexion

An additional observation from FIG. 53 is that when testing at a bias of0 degrees, the elastic stiffness increases from trial 1 to trial 10. Asstated previously, trials 1 and 2 have a maximum frequency component of4 Hz, trials 3 and 4 have a maximum frequency component of 5 Hz, trials5 and 6 have a maximum frequency component of 6 Hz, and trials 7 and 8have a maximum frequency component of 8 Hz. One hypothesis that canexplain this observation is that because velocities of the perturbationsincrease with trial number, the primary afferents of the muscle spindlesense these greater rates of muscle stretch in turn activating thestretch reflex causing muscles to contract and become stiffer. Thus, thehypothesis says that the combined active stiffness of the muscles actingon the wrist may be increasing with trial number. An observed increasein EMG activity at the greater trial numbers supports this hypothesis.This trend, however, is not shown at the bias of 30. This is likelybecause extensor muscles are already active in all 10 trials as they arestretched to a length further from the neutral position.

Many subjects with spasticity cannot be tested at a neutral biasposition of 0 degrees as they have difficulty extending their wristjoint to that position. The resting configuration of the hand of someonewith spasticity is in the hyper-flexed position due to hyperactivity ofthe flexor muscle group. Four adult patients with hemiplegic spasticitywere tested, AM, JS, SH, and BY, 3 males and 1 female. A hemiplegicpatient is one where one side of his or her body is affected byspasticity whereas the contra-lateral side is practically unaffected orless affected. Because the contra-lateral arm is unaffected or lessaffected, it can serve as a control to the spastic arm when testing apatient with hemiplegic spasticity. JS is a patient with Parkinson'sdisease, SH suffers from traumatic brain injury due to a motorcycleaccident, and BY is a patient who severely lost function of his lefthand after a stroke. AM's affected arm was tested at a bias of 30degrees flexion. For subjects JS, SH, and BY the affected arms andunaffected contra-lateral arms were tested, both at a bias of 30degrees. The results of these experiments are shown in FIGS. 54 and 55.Unlike the figures showing results of the control subjects, all plotsare not on the same scale due to the large range in values. The whitebar graphs show results for the unaffected arm and the grey bar graphsshow results for the contra-lateral affected arm. Both arms were testedat a bias of 30 degrees flexion.

Looking at the average of {tilde over (τ)}_(off), {tilde over (K)}_(H)and {tilde over (B)}_(H) over all ten trials for each experiment denotedby the asterisk over the eleventh bar labeled trial 11, the followingconclusions can be made. In the experiment with the unaffected hand atthe bias of 30 degrees flexion, the average {tilde over (τ)}_(off)across all ten trials are comparable to what is seen in control subjectswhen tested at the same bias. As in the control subjects, the torqueoffset values are positive indicating that the unaffected hand prefersto be at the neutral 0 degree bias position as in the control subjects.However, when testing the affected, spastic side at the 30-degree bias,the {tilde over (τ)}_(off) values are negative in sign and greater inmagnitude. This indicates that the preferred resting position of thesubject's spastic wrist joint is at a flexed position even greater than30 degrees, as the hand is pushing towards an even further flexedposition. This phenomenon is what is called hyper-flexion. Thisobservation of increase in magnitude and change in sign in the average{tilde over (τ)}_(off) across the ten trials when testing the affectedhand versus the unaffected hand is significant (p=0.0177). All 3subjects tested on both arms have a higher average stiffness (mean of{tilde over (K)}_(H) across the ten trials) in their spastic arm ascompared to their contra-lateral arm (p=0.0743). When comparing theaverage stiffness between contralateral arms, in subjects JS, and SH,the difference is large. Also, in general, the average viscosity, {tildeover (B)}_(H), across the ten trials is greater in the affected handwhen compared to the unaffected hand, but the difference seen inviscosity values is not as great as the difference seen in the stiffnessvalues (p=0.199). Here, the difference between the affected spastic armis compared with the unaffected contra lateral arm within the spasticgroup.

The average and standard deviation (denoted by the height of the errorbars) of {tilde over (τ)}_(off), {tilde over (K)}_(H) and {tilde over(B)}_(H) values across all control adults for each trail are shown inFIG. 56. The eleventh bar shows the overall average and standarddeviation across all control adults, across all ten trials.

Table 5 below shows the average and standard deviation of {tilde over(τ)}_(off), {tilde over (K)}_(H) and {tilde over (B)}_(H) values acrossall ten FIG. 56. Table 6 below shows the range of average {tilde over(τ)}_(off), {tilde over (K)}_(H) and {tilde over (B)}_(H) values acrossall ten trials in the group of adults with hemiplegic spasticity.

TABLE 5 Hand tested {tilde over (τ)}_(off) [N * m]: {tilde over (K)}_(H){tilde over (B)}_(H) (both at bias of 30 Torque offset [N * m/rad]: [N *m * s/rad]: degrees flexion) towards extension Stiffness ViscosityUnaffected 0.213 ± 0.0877 2.062 ± 0.820 0.0709 ± 0.0318 Affected −0.593± 0.471    3.858 ± 1.648 0.0992 ± 0.0447

TABLE 6 {tilde over (τ)}_(off) [N * m]: Hand tested Torque offset {tildeover (K)}_(H) (both at bias of 30 towards [N * m/rad]: {tilde over(B)}_(H) [N * m * s/rad]: degrees flexion) extension Stiffness ViscosityUnaffected 0.144 to 0.281 1.220 to 2.857 0.0343 to 0.0916 Affected−1.1008 to 2.606 to 6.284 0.0694 to 0.166  −0.161

Tables 5 and 6 show that at the bias position of 30 degrees flexion, theaverage parameter values of the unaffected or less affected arms of thehemiplegic population are still somewhat higher than when compared tothe control population. This demonstrates that the tone of the lessaffected side of the patient group in the is still greater on averagethan what would be seen in the control group, though this difference isof course greater when comparing the severely affected side of thehemiplegic population with the control population. This is a fairercomparison since subjects in the control group have no knownneurological disorders. When comparing the hemiplegic subjects to thecontrols, the observation of the negation of sign and increase inmagnitude of the average torque offset values, {tilde over (τ)}_(off)across the ten trials shown by the 11^(th) bars are extremelysignificant (p=3.0*10⁻¹⁰). The increase in average stiffness, {tildeover (K)}_(H) and viscosity, {tilde over (B)}_(H) values across the tentrials when compared across populations are also extremely significant(p=1.5*10⁻⁸ and p=6.5*10⁻⁵ respectively).

On average, the stiffness and viscosity values across all trials of thecontrol group were 1.305 N*m/rad and 0.0512 N*m*s/rad respectively. Thecorresponding stiffness and viscosity values of the severely affectedside of the hemiplegic population were 3.858 N*m/rad and 0.0992N*m*s/rad respectively. These elevated values when compared to those ofthe control group are extremely significant (p=1.5*10⁻⁸ for stiffnessand p=6.5*10⁻⁵ for viscosity).

While the device and method have been described in detail for use inquantifying muscle tone, particularly in the wrist, in a spasticpatient, the invention is not so limited. Rather, the device and methodare useful on both the upper and lower extremities including, forexample, the ankle. Thus, for example, the device could be modified tosuit the leg and ankle accordingly. Further, the method and device canbe used to measure muscle tone in the non-spastic patient. For example,the method and device can be used to measure rigidity in a patient.Further, the method and device could be used in analyzing a patient'smuscle tone in line with, for example, physical therapy that the patientis undergoing to regain control of the muscles in the legs after aspinal accident.

Although a preferred embodiment of the invention has been describedusing specific terms, such description is for illustrative purposesonly, and it is to be understood that changes and variations may be madewithout departing from the spirit or scope of the following claims.

What is claimed is:
 1. A method for quantifying muscle tone in a patientcomprising the steps of: moving the patient's wrist in a non-sinusoidaland non ramp trajectory and determining the stiffness, viscosity andinertial parameters using the following relationship${\tau_{s}(t)} = {{K_{H}{\theta (t)}} + {B_{H}{\overset{.}{\theta}(t)}} + {J_{T}{\overset{¨}{\theta}(t)}} + \tau_{off}}$

where τ_(s) is the total torque, τ_(off) is the offset torque, K_(H) andB_(H) are the angular stiffness and viscosity of the combined flexor andextensor muscle groups that act on the joint, J_(T) is the combinedinertia of the oscillating appendage, θ is the angular displacement ofthe system, and$\overset{.}{\theta}\quad {and}\quad \overset{¨}{\theta}$

 are the velocity and acceleration.
 2. A method of quantifying muscletone in a patient comprising the steps of: driving the patient's wristis in arbitrary trajectory that is neither sinusoidal nor ramp;controlling the trajectory θ(t); measuring the torque response τ(t); anddetermining the stiffness K_(H), viscosity B_(H), and inertial J_(T)parameters using the relationship${\tau_{s}(t)} = {{K_{H}{\theta (t)}} + {B_{H}{\overset{.}{\theta}(t)}} + {J_{T}{\overset{.}{\theta}(t)}} + \tau_{off}}$

wherein where τ_(s) is the total torque, τ_(off) is the offset torque,K_(H) and B_(H) are the angular stiffness and viscosity of the combinedflexor and extensor muscle groups that act on the joint, J_(T) is thecombined inertia of the oscillating appendage, θ is the angulardisplacement of the system, and$\overset{.}{\theta}\quad {and}\quad \overset{¨}{\theta}$

 are the velocity and acceleration.
 3. A device for quantifying muscletone in a patient comprising: an electromechanical system that drivesthe wrist of a patient with an arbitrary trajectory that is neithersinusoidal nor ramp, wherein the system includes: a housing having a topsurface; a hand support on the top surface of the housing for holdingthe hand in a desired position; an arm support on the top surface of thehousing for holding the arm in a desired position relative to the hand;a motor coupled to the hand support such that the motor moves back andfourth tracking the desired trajectory and wherein the motor causes thehand to move in the desired trajectory; and an automated feedbackcontroller for varying the torque, speed and direction of the motor. 4.The device of claim 3, wherein a discrete set of finite data pointsdefining the trajectory is input into the controller, and wherein themotor tracks these data points.
 5. The device of claim 3, wherein thehand securing mechanism comprises at least a pair of vertical barsextending upwards from the top surface of the housing and being spacedapart so the hand can be placed between the pair of vertical bars. 6.The device of claim 5, wherein the vertical bars have padding or acushion on their outer surface.
 7. The device of claim 6, wherein thepadding or cushion is removable and varies in thickness for adapting thedevice to various sized hands.
 8. The device of claim 6, wherein thepadding or cushion is compressible and expandable to accommodate handsof varying sizes.
 9. The device of claim 5, wherein the vertical barsare slidably attached to the housing so that the vertical bars aremovable towards each other and away from each other and locked into adesired position for accommodating varying sized hands.
 10. The deviceof claim 3, further comprising an arm securing mechanism to secure thearm in place.
 11. The device of claim 3 further comprising an opticalencoder for converting motion of the device and wrist into a series ofdigital pulses.
 12. The device of claim 3 further comprising anamplifier for amplifying the command voltage from the controller intoeither a voltage or a current to the motor.